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ABSTRACT 


The  objective  of  this  study  is  the  development  of  a 
procedure  for  determining  the  large  dynamic  response  of 
structural  systems  consisting  of  beam  and  rectangular 
plate  elements.  The  analysis  takes  into  account  both 
geometrical  and  material  nonlinearities.  The  general 
approach  to  the  problem  is  based  on  the  finite  element 
method  and  the  use  of  displacement  interpolation  functions. 
The  energy  expressions  for  both  the  beam  and  plate 
elements  are  obtained  and  are  used  to  generate  the  stiffness 
and  mass  matrices  of  the  elements.  Also,  the  geometric 
stiffness  matrices  are  derived  which  account  for  the 
effect  of  geometric  nonlinearities.  Plastic  deformations 
are  taken  into  account  by  means  of  an  incremental  theory 
of  plasticity  coupled  with  the  concept  of  initial  strain. 

A  computer  program  is  developed  for  the  analysis  of 
structures  which  consist  of  beam  and  rectangular  plate 
elements.  This  program  may  be  used  to  perform  an  elastic 
static,  elastic  dynamic,  or  plastic  dynamic  analysis  with 
or  without  the  inclusion  of  geometric  nonlinearity  effects. 
It  can  also  be  used  to  perform  a  free  vibration  analysis 
of  a  structure,  leading  to  the  natural  frequencies  and  modes 
of  vibration.  The  computer  program  is  used  in  this  study 
to  obtain  the  solution  to  several  example  structures 


subjected  to  static  or  dynamic  loads.  Results  indicate  that 
both  geometric  and  material  nonlinearities  have  an  important 
effect  in  the  deformation  of  structures  composed  of  beam 
and  plate  elements. 
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I .  INTRODUCTION 


1.1  General  and  Review  of  Literature 

Large  deformations  of  a  structural  system  give  rise  to  a 
nonlinear  phenomenon.  This  leads  to  nonlinear  equations  which 
immediately  render  classical  methods  of  analysis  inapplicable. 
Among  many  approaches  to  the  solution  of  the  problem,  the 
finite  element  method  has  been  the  one  most  commonly  used  in 
recent  years.  The  structure  is  first  idealized  as  an  assembly 
of  discrete  structural  elements,  and  some  assumptions  on  the 
displacement  distribution  with  respect  to  these  elements  are 
then  made.  The  complete  formulation  is  obtained  by  combining 
these  individual  approximate  displacement  distributions  in  a 
manner  which  satisfies  the  force-equilibrium  and  displacement 
compatibility  relations  at  the  junctions  of  these  elements. 
Numerical  methods  of  analysis  may  then  be  used  to  solve  the 
equations.  In  such  a  formulation  of  the  problem,  matrix  algebra 
is  appropriately  utilized  and  conveniently  lends  itself  to  solu¬ 
tion  by  means  of  a  digital  computer. 

In  dealing  with  large  deformations  of  structures  two  types 
of  nonlinearities  have  to  be  taken  into  account,  i.e.,  geometrical 
and  material  nonlinearities.  The  concept  of  using  the  finite 
element  approach  in  solving  geometrically  nonlinear  structures 
was  initially  advanced  by  Turner  et  al  (1) .  Martin  (2)  reviewed 
the  work  on  geometrically  nonlinear  problems  and  presented  geo¬ 
metric  stiffness  matrices  for  a  number  of  structural  elements, 


while  Oden  (3)  gave  special  attention  to  incremental  formulations 
of  these  problems.  Oden  also  studied  the  problem  from  dual  view¬ 
points:  first,  the  approach  in  which  nonlinear  stiffness  relations 

with  complete  generality  were  generated  and  solved  by  the  Newton- 
Raphson  method.  Secondly,  the  approach  in  which  the  structural 
problem  was  viewed  as  one  of  minimizing  a  scalar-valued  function 
defining  the  total  potential  energy  of  the  system  of  n-variables. 

The  behavior  of  geometrically  nonlinear  frames  was  the 
subject  of  several  investigations,  such  as  (4-7) .  The  problem 
of  large  deformation  of  plate  structures  was  previously  investi¬ 
gated  by  Timoshenko  (8) ,  who  considered  the  interaction  of 
bending  and  in-plane  (membrane)  action.  Greene  (9) ,  and  Kapur  (10) 
derived  separately  the  stiffness  properties  for  a  rectangular 
plate  for  various  static  loading  conditions.  Murray  and  Wilson  (11) 
studied  the  large  deflection  of  plates  by  using  triangular  elements 
in  which  in-plane  as  well  as  out-of-plane  displacements  were 
permissible.  Brebbia  and  Connor  (12)  investigated  the  same  prob¬ 
lem  by  initially  setting  up  the  stiffness  matrix  for  the  structure 
and  then  solving  the  system  of  force-displacement  equations  using 
a  limited  number  of  load  steps,  and  followed  by  corrections  based 
on  the  Newton-Raphson  method. 

Zienkiewicz  and  Cheung  (13)  analyzed  a  plate  with  edge  beams 
subjected  to  static  loading  of  relatively  small  magnitude,  and 
compared  the  finite  element  solution  to  the  classical  solution 
given  by  Timoshenko  (8) .  The  elastic-plastic  analysis  of  plates 
with  edge  beams  was  considered  by  McNeice  (14) .  Weaver  et.al.  (15) 


considered  a  rectangular  plate  with  in-plane  and  torsional  action. 
Timoshenko  (16)  studied  the  problem  of  the  small  deformation  of 
a  rectangular  plate  subject  to  a  concentrated  lateral  load  at 
center  and  uniformly  compressed  in  the  in-plane  directions. 

In  recent  years,  there  have  been  other  significant  contri¬ 
butions,  for  example  (17,18).  Furthermore,  the  problem  of 
material  nonlinearity  as  related  to  structural  systems  has 
received  increased  attention  by  several  investigators,  such  as 
(19-21).  In  most  of  these  investigations,  plastic  behavior 
beyond  the  elastic  range  has  been  determined  by  the  use  of  the 
Prandtl-Reuss  equations  (22)  of  the  classical  plasticity  theory. 
Concurrently,  two  distinct  methods  of  incorporating  plastic  behavior 
into  a  finite  element  analysis  have  been  developed,  and  are  refer¬ 
red  to  in  the  literature  as  the  "initial  strain"  and  the  "tangent 
modulus"  methods  (20) .  Such  techniques  are  very  suitable  in 
dealing  with  problems  that  involve  both  geometric  and  material 
nonlinearities . 

Very  few  solutions  to  structural  dynamics  problems  have 
been  reported  in  the  literature  which  aim  at  determining  the 
complete  response  of  a  structural  system,  including  the  elastic 
and  inelastic  deformations,  as  well  as  the  effect  of  geometric 
nonlinearities.  In  this  connection  Farhoomand,  Iverson,  and 
Wen  (23)  developed  a  method  of  dynamic  analysis  of  space  frames 
considering  the  combined  effect  of  material  and  geometric  non- 
linearities.  In  their  analysis,  the  mass  properties  of  the 
structural  components  were  lumped  at  the  joints  and  the  material 
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was  assumed  to  be  elastic-perfectly  plastic,  with  inelastic 
deformations  taking  place  at  member  ends  only.  Likewise, 

Toridis  et  al  (24,25)  presented  an  incremental  finite  element 
procedure  for  the  dynamic  analysis  of  framed  structures  with  both 
geometric  and  material  nonlinearities.  Their  solution  was  based 
on  the  concept  of  initial  strain  and  the  use  of  the  geometric 
stiffness  matrix  for  a  beam’ element.  Strain  hardening  and 
hysteretic  effects  were  also  taken  into  account. 

More  recently,  McNamara  and  Marcal  (26)  developed  a  numerical 
procedure  aimed  at  the  use  of  the  finite  element  method  to  analyze 
large  elastic-plastic  deformations  of  structures  under  a  variety 
of  dynamic  loadings.  In  their  procedure,  the  basic  equations 
for  nonlinear  finite  element  analysis  were  linearized  based  on 
an  incremental  process  and  the  addition  of  two  error  terms  to 
account  for  discretization  errors.  It  was  shown  that  the 
corrected  incremental  equations ,  in  the  form  used  in  their 
analysis  give  stable  and  accurate  solutions  even  with 
relative  large  time  increments. 

1.2  Object  and  Scope 

The  objective  of  the  present  investigation  has  been  the 
development  of  a  procedure  for  determining  the  large  dynamic 
response  of  structural  systems  consisting  of  beam  and  rectangular 
plate  elements.  The  analysis  takes  into  account  both  geometrical 
and  material  nonlinearities.  The  general  approach  to  the  problem 
is  based  on  the  finite  element  method  and  the  incremental  plasticity 
theory  coupled  with  the  use  of  the  geometric  stiffness  matrices  of 
the  elements  to  account  for  significant  changes  in  the  geometry  of 


the  structure.  From  the  practical  viewpoint,  one  of  the  additional 
objectives  of  this  study  has  been  the  development  of  a  general 
purpose  computer  program  for  the  analysis  of  the  type  of  structures 
under  consideration.  The  computer  program  is  designed  with  suf¬ 
ficient  flexibility  for  possible  use  with  various  classes  of 
structural  elements.  Presently,  the  classes  of  elements  that 
have  been  incorporated  into  the  program  consist  of  two  types  of 
beam  elements  corresponding  to  plane  or  three-dimensional  action, 
and  five  types  of  thin  rectangular  plate  elements  corresponding 
to  membrane  and/or  bending  and  twisting  action.  The  beam  and 
plate  elements  may  be  used  alone  or  in  combination  (beam-plate 
assemblage)  to  perform  a  static,  elastic  dynamic  or  free  vibration 
analysis  of  a  structural  system.  Likewise,  a  plastic  dynamic 
analysis  may  be  performed,  but  the  type  of  plate  elements  that 
can  be  used  for  this  purpose  is  presently  restricted  to  rectangular 
elements  with  membrane  (in-plane)  action  only. 
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II.  BRIEF  REVIEW  OF  BASIC  CONCEPTS 


The  use  of  the  finite  element  method  in  solving  structural 
problems  is  well  known.  For  comprehensiveness,  related  topics 
such  as  Castigliano1 s  theorem,  the  strain  energy  expression, 
and  the  stress-strain  relationship  for  the  problem  under  con¬ 
sideration  are  briefly  discussed.  The  concept  of  the  direct 
stiffness  method  is  next  outlined  and  is  followed  by  the  formula¬ 
tion  of  the  dynamical  expressions  governing  the  behavior  of  the 
system.  Temperature  is  assumed  to  be  constant  throughout  all 
the  derivations. 

2.1  Brief  Review  of  Basic  Concepts 

2.1.1  Strains  and  Displacements 

The  structure  under  study  is  first  partitioned  into 
a  set  of  idealized  elements.  For  each  idealized  element,  the 
interior  displacements  [u]  are  assumed  to  be  expressible  in  terms 
of  the  nodal  displacements  [U]  by  the  matrix  equation 

[u]  =  [a]  [U]  (1) 

where  [a]  =  [a(x,y,z)]  is  a  rectangular  matrix  and  a  function 
of  the  position  coordinates.  The  matrix  [a]  is  usually  determined 
by  assuming  the  form  of  the  displacement  distribution.  This 
displacement  distribution  must  be  continuous  and  should  preferably 
satisfy  compatibility  of  deflections  and  slopes  on  the  boundaries 
of  the  elements.  The  satisfaction  of  the  stress-equilibrium 
equations  is  also  desirable.  It  is  thus  clear  that  in  general, 
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only  approximate  expressions  for  the  matrix  [a]  are  expected 
since  the  latter  depends  on  the  complete  displacement  distri¬ 
bution  of  the  entire  structure. 

Since  the  total  strains  [e]  are  functions  of  the  dis¬ 
placements  [u] ,  which  in  turn  are  functions  of  [U] ,  the  follow¬ 
ing  matrix  equation  is  valid: 

[e]  =  [b]  [U]  (2) 


where  [b]  is  again  a  rectangular  matrix. 

2.1.2  Castigliano 1 s  Theorem,  Strain  Energy  Expression 
Stress-Strain  Relationship 

The  well-known  Castigliano ' s  Theorem  (27)  states  that 

if  a  structure  is  subjected  to  a  system  of  external  forces 

S-,  ,  S„,  ...,  S  and  if  only  one  virtual  displacement  6u  is 
J.  ^  n  y 

applied  in  the  direction  of  the  load  S^,  and  if  V  represents  the 
strain  energy  of  the  system,  then 


(3) 


If  Uy  is  identified  as  a  nodal  displacement  and  if  the  index  is 
varied  from  y  =  1  to  y  =  n,  then  the  following  matrix  equation 
may  be  obtained  from  this  set  of  n  equations 


[S] 


9  V 

TTuT 


(4) 


The  strain  energy,  V,  shown  in  Eq.  (4) ,  may  be  expressed 
in  matrix  notation  as  (dv  represents  an  infinitesimal  volume  element) 

1  T 

V  =  Tfv  [E]  ta]  dv*  (5) 
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The  stresses  [a]  and  strains  [e]  shown  in  Eq.  (5)  are 
related  by  the  following  equation: 

[a]  =  [D]  [el  (6) 

where  [D]  is  a  rectangular  matrix,  referred  to  as  the  "elasticity 
matrix." 

2.2  Stiffness  and  Inertial  Properties  of  Elements 

The  concept  underlying  the  direct  stiffness  method  is  the 
determination  of  the  displacements  of  an  idealized  structure 
under  some  specified  external  loading  by  investigating  the 
stiffness  properties  of  individual  structural  elements  which  can 
be  used  to  represent  the  idealized  structure.  The  displacements 
will  be  considered  as  the  unknowns  throughout  the  analysis.  The 
stiffness  properties  of  an  individual  element  may  be  obtained 
by  substituting  Eqs.  (2) ,  (6)  into  Eq.  (5) ,  so  that 

T 

V  =  j  fv  [e]  [a]  dv  =  \  /v  [U]T  ( [b] T  [D]  [b] )  [U]dv.  (7). 

Differentiating  the  above  equation  with  respect  to  the  ith  nodal 
displacement  and  substituting  the  result  into  Eq.  (4) ,  the  follow¬ 
ing  equation  is  obtained 

[S]  =  =  (/v  (  [b] T  [D]  [b]dv)  [U]  .  (8) 

Defining 

[K]  =  fv  ( [b]T  [D]  [b]dv  (9) 

where  [K]  represents  the  elastic  stiffness  matrix  of  the  element, 
Eq.  (8)  now  becomes 
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ts]  =  [K][U] 


(9a) 


which  is  the  force-deformation  relationship  of  the  element  under 
consideration.  It  can  be  observed  from  Eqs.  (8) -(9a)  that  the 
stiffness  matrix  components  may  be  represented  as 


2 

jk  3U.3U, 

J  j  k 


(10) 


It  should  also  be  observed  that  in  view  of  Eqs.  (7)  and  (9) ,  the 
strain  energy  expression  V  may  be  represented  as 


V  =  \  [U]T[K]  [U] 


(11) 


If  only  one  element  is  used  to  idealize  the  structure, 

Eq.  (9a)  is  the  representative  equation.  However,  in  the  general 
case  of  a  structure  consisting  of  several  elements,  the  force- 
deformation  relation  of  all  the  elements,  as  given  in  Eq.  (9a), 
may  be  appropriately  assembled  to  form  the  representative  equa¬ 
tion  of  the  entire  system. 

In  considering  the  dynamic  response  of  a  structure,  the 
inertial  properties  of  the  system  must  first  be  determined.  As 
shown  in  Refs.  (25)  and  (27)  ,  the  mass  matrices  of  the  individual 
elements  of  the  structure  may  be  computed  through  a  process 
similar  to  the  generation  of  the  stiffness  matrices.  Again,  an 
integration  over  the  volume  of  the  element  is  needed,  so  that, 
letting  [M]  represent  the  element  mass  matrix 

[M]  =  /v  p  [a]  T  [a]  dv  (12) 


where,  p,  represents  the  mass  density  of  the  material. 
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As  in  the  case  of  the  stiffness  matrices,  the  mass  matrices 
of  the  individual  elements  may  be  assembled  appropriately  to  form 
the  generalized  mass  matrix  for  the  entire  structure. 

2.3  Effect  of  Geometric  and  Material  Nonlinearities 

As  will  be  shown  later  in  this  report,  the  effect  of 
significant  changes  in  the  geometry  of  the  structure  can  be 
accounted  for  by  the  inclusion  of  higher  order  nonlinear  terms 
in  the  strain-displacement  relations.  This,  in  turn,  intro¬ 
duces  an  additional  stiffness  matrix  K_,,  referred  to  as  the 

(j 

"geometric  stiffness  matrix."  In  such  a  case,  it  has  already 
been  demonstrated  in  Ref.  (25)  that  for  a  framed  structure 
consisting  of  beam  elements,  the  strain  energy  expression  given 
by  Eq.  (12)  takes  the  form  of 

V  =  i  [U]T([K]  +  [Kg]  )  [U]  (13) 

where,  as  before,  [K]  represents  the  elastic  (first  order) 
stiffness  matrix  of  the  element  defined  by  Eq.  (9) ,  and  [Kg] 
refers  to  the  geometric  stiffness  matrix. 

When  inelastic  deformations  take  place,  it  is  assumed 
(24)  that  the  total  strain  vector  [ e ]  may  be  divided  into  an 
elastic  part  [e  ]  and  a  plastic  part  [er] .  In  turn,  the 
plastic  strain  vector  may  be  treated  as  an  initial  strain  vector 
[e°] ,  so  that 

[eG]  =  [e]  -  [e°]  (14) 
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In  a  similar  manner,  it  is  assumed  that  the  nodal  displacements 
of  the  element  may  be  separated  into  an  elastic  part  [Ue]  and 
plastic  or  initial  part  [U°] .  Furthermore,  it  is  assumed  that 
the  interpolation  or  shape  functions  chosen  to  represent  the 
elastic  displacement  distributions  within  the  elements  may  also 
be  used  to  represent  the  plastic  displacement  distribution  in  the 
elements.  Following  this  approach,  it  is  shown  in  Refs.  (24) 
and  (25)  that  the  effect  of  the  plastic  actions  occurring  in  the 
elements  may  be  represented  as  equivalent  forces  acting  at  the 
nodal  points  of  the  elements.  In  this  manner,  the  modified 
equilibrium  equations  can  be  easily  formulated. 

2.4  Transformation  Matrices 

So  far,  the  discussion  has  centered  on  the  generation  of 
the  stiffness  and  mass  properties  of  an  individual  element 
referred  to  the  local  coordinate  system  of  the  element.  In 
order  to  determine  the  corresponding  properties  of  the  entire 
structure,  a  global  coordinate  system  must  be  established.  The 
displacements  and  corresponding  forces  of  the  elements  are  then 
transformed  to  the  global  coordinate  system,  prior  to  the  appli¬ 
cation  of  the  direct  stiffness  method. 

Let  [U]  and  [S]  be  the  element  displacements  and  forces 
with  respect  to  the  local  coordinate  system,  and  [U]  and  [S]  be 
the  same  element  displacements  and  forces  with  respect  to  the 
global  coordinate  system,  respectively.  Let  [X]  be  the  trans¬ 
formation  matrix  between  these  two  coordinate  systems.  Then, 

[U]  =  [X]  [U].  (15) 
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Additionally,  by  using  the  principle  of  virtual  work,  the  follow¬ 
ing  expression  may  be  written 

[S]  =  [X]  [s]  .  (16) 

Substitution  of  these  two  equations  into  Eq.  (10)  yields 

[S]  =  [K] [U]  (17) 

where 

[K]  =  [X]T[K][X]  (18) 

Eq.  (17)  is  analogous  to  Eq.  (10)  except  that  the  stiffness 
matrix  of  anyone  element  is  now  expressed  with  reference  to  the 
global  coordinate  system. 

The  element  mass  matrices  may  be  transformed  in  a  similar 
manner.  Letting  [M]  represent  the  element  mass  matrix  with 
respect  to  the  global  coordinate  system 

[M]  =  [A]T[M]  [X]  (19) 

2.5  The  Basic  Dynamical  Equation 

The  dynamical  equation  of  the  system  is  derived  based  on 
Hamilton's  Principle  and  the  calculus  of  variations.  The 
Hamiltonian  function  is  formed  by  considering  the  total  kinetic 
and  potential  energies  of  the  system.  The  variation  of  the 
integral  of  the  Hamiltonian  function  between  two  discrete  times 
t^  and  t2  is  then  set  equal  to  zero  in  accordance  with  Hamilton's 
Principle.  The  basic  dynamical  equation  for  the  structure  under 
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consideration  is  then  derived  by  following  through  some  rather 
tedious  mathematical  operations,  as  reported  in  Refs.  (24)  and 
(25) .  Since  details  of  the  various  steps  involved  in  this 
derivation  are  given  in  the  above  references,  they  will  not  be 
repeated  here.  The  resulting  expression  is  in  the  form  given 
below: 


[m]  [q]  +  (  [k]  +  [kG]  )  [q]  =  [F]  +  [F°] 


(20) 


where 


[m] 

[q],[q] 

[k] 

[kG] 

[F] 

[F°] 


assembled  (generalized)  mass  matrix. 

generalized  displacement  and  acceleration  vectors, 
respectively. 

assembled  (generalized)  elastic  stiffness  matrix, 
assembled  (generalized)  geometric  stiffness  matrix, 
generalized  nodal  force  vector. 

( [k]  +  [kg])  [q°]  =  equivalent  generalized  nodal  force 

vector  due  to  plastic  strains. 


Eq.  (20)  represents  the  basic  equation  governing  the  dynamic 
behavior  of  the  entire  structure.  In  considering  the  effect  of 
significant  changes  in  the  geometry  of  the  structure  due  to  large 
deformations  (second  order  analysis) ,  both  the  stiffness  and  mass 
matrices  have  to  be  recomputed  at  various  time  intervals,  based 
on  the  current  deformed  configuration  of  the  structure.  On  the 
other  hand,  eliminating  the  matrix  [k  ]  from  Eq.  (20)  leads  to 
the  basic  equation  for  the  small  deformation  (first  order)  analysis 
of  the  structure.  Thus,  it  is  clear  that  the  particular  expressions 
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corresponding  to  the  first  and  second  order  static,  elastic 
dynamic  and  plastic  dynamic  analysis  of  the  structure  under 
consideration  may  be  obtained  easily  from  Eq.  (20)  by  con¬ 
sidering  only  those  terms  in  the  equation  that  are  appropriate 
in  each  case. 


III.  THE  BEAM  ELEMENT 


The  fundamental  relations  outlined  in  the  previous  chapter 
form  the  basis  of  the  present  study  aimed  at  analyzing  the 
behavior  of  a  structural  system  consisting  of  straight  beam  and 
rectangular  thin  plate  elements.  In  forming  these  expressions, 
the  properties  of  the  individual  elements  of  the  structure  must 
first  be  derived.  In  this  chapter  the  force-deformation  and 
inertial  properties  of  a  three-dimensional  beam  member  are  dis¬ 
cussed  and  the  stiffness  and  mass  matrices  are  presented.  The 
behavior  of  a  rectangular  thin  plate  element  is  studied  in  the 
next  chapter. 

3.1  Strain  Energy  Expression  for  a  Beam  Element 

Consider  a  typical  beam  element  with  an  arbitrary  orienta¬ 
tion  in  space,  as  shown  in  Fig.  1.  A  set  of  axes  X,  Y,  Z  are 
chosen  to  represent  a  global  coordinate  system  of  reference. 

The  element  local  coordinate  axes  are  x,  y,  z  with  the  x-axis 
extending  along  the  axis  of  the  element,  and  the  y,  z-axes 
pointing  in  the  principal  directions  of  the  cross  section.  The 
set  of  element  nodal  (end)  displacements  and  forces  are  designated 
by 


[U]  =  [U1,U2,...,U12]  and  [S]  =  [S]L  ,S2 , .  .  .  ,S121  , 

respectively. 

In  order  to  define  the  element  displacements  and  forces  with 
respect  to  the  global  system  of  axes,  both  during  the  initial  and 
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deformed  configurations  of  the  elements,  a  transformation  matrix 
[X]  is  utilized  (see  Eq.  (15)).  Letting  x,y,z  represent  the 
current  coordinates  of  an  element  nodal  point  (the  original  co¬ 
ordinates  of  the  nodal  point  plus  the  current  translational 
displacements  of  the  point) ,  the  transformation  relation  is  expressed 

X 

[A]  Y  (21) 

Z 

in  which  [X]  is  represented  as  a  diagonal  matrix,  so  that  in 
partitioned  form 


[X]  =  diag  ([T^] , [T2] , [T3] , [T4])  (22) 

where 

[T1]  =  [T2]  =  [T3]  =  [T4]  = 

In  the  above  m^  and  n^  are  the  instantaneous  direction 

cosines  of  the  ith  local  coordinate  axis  with  reference  to  the 

X,  Y  and  Z  global  axes,  respectively. 

Let  [u]  represent  the  vector  of  element  displacements  at  any 

point  in  the  interior  of  the  element  with  components  ux,  u^  and 

u  along  the  x,  y  and  z-directions ,  respectively.  Making  use 
z 

of  simple  beam  theory,  [u]  can  be  expressed  in  terms  of  the  nodal 
displacement  vector  [U]  as 
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[u] 


[a]  [U] 


(24) 


where,  the  transpose  of  the  matrix  [a] ,  containing  the  inter¬ 
polation  (shape)  functions,  has  the  form  of 
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When  large  deformations  take  place,  the  normal  strain  com¬ 


ponent  e  in  the  local  x-direction  contains  certain  non-linear 

XX 


terms  which  cannot  be  neglected.  In  terms  of  the  element  dis¬ 


placements,  e  may  be  represented  as 

XX 


3ux  32u 
exx  =  “at  -y 


.  3u  32u  i 

1  +  1  (_Z)2  -z _ 2  1  __z 

2  +  2  3x  „  2  2  1  3x' 


3u_  2 


(26) 


3x  3x 

In  this  expression  the  first,  second  and  fourth  terms  represent 
the  contributions  of  the  axial  and  bending  forces  in  the  usual 
small  deflection  theory.  The  third  and  fifth  terms  are  approxi¬ 
mations  to  the  contributions  due  to  rotations  resulting  from 
bending  around  the  z  and  y-axes,  respectively. 

In  order  to  obtain  initial  estimates  of  the  displacements 
within  each  increment  of  loading,  the  elastic  relations  based 
on  Hooke's  Law  may  be  utilized  at  the  beginning  of  each  time  step. 
The  strain  energy  in  a  beam  element  may  then  be  expressed  as 


,T  Ef  2  ,  ,  G  ,  2  , 

V  =  /,T  e  dv  +  7r  /  e  dv 

2  v  xx  2  yz 


where  E  is  the  modulus  of  elasticity. 


(27) 

Combining  Eqs.  (26)  and 


(27) 


„  3u 

V  =  —  /  [ - —  —  y. 

2  v  1  3x  y 


3  2u 


u 


i  2 

1  +  l  (_X)  - 

22  (  dx1 


ri  3U 

+  2  [  (t£  + 


32u  i  2  2 

:— |  +  i  (-sf)2i 

3x 
3u 


3y 


)]  dv 


(28) 


Each  individual  term  in  the  above  brackets  can  be  evaluated 
based  on  Eqs.  (24)  and  (25).  Substituting  the  resulting  relations 
into  the  expression  for  V,  neglecting  higher  order  terms  (such  as. 
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4  i  9uz  4 

,  (~g^)  )  and  performing  the  indicated  integration, 

Eq.  (28)  can  then  be  rewritten  in  the  form  given  previously  as 
Eq.  (13). 

3.2  Stiffness  and  Mass  Matrices  for  a  Beam  Element 

Referring  to  Eq.  (13)  ,  for  a  beam  element,  the  elastic 
matrix  [K]  and  geometric  stiffness  matrix  [K„]  have  been 
derived  in  Refs.  (25)  and  (27)  ,  and  are  in  the  form  given 
below. 
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where,  in  Eqs.  (29)  and  (30),  GI  is  the  torsional  stiffness 
of  the  beam  cross  section  and  I  and  I  refer  to  the  moments 

y  z 

of  inertia  around  the  y  and  z  axes,  respectively.  A  denotes  the 
cross  sectional  area  of  the  beam  element,  H  is  its  length  and 
G  represents  the  shear  modulus.  The  axial  force  N  is  defined  as 


f\,  FA 

N  =  constant  =  —j  (U^  -  U^) 


in  which,  according  to  the  previous  definition,  and  refer 
to  the  nodal  displacements  in  the  axial  direction,  corresponding 
to  the  left  (first)  and  right  (second)  end  of  the  beam  element, 
respectively. 

Based  on  Eq.  (24)  the  velocity  and  acceleration  distribution 
within  the  element  may  be  approximated  as 
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[u] 

=  [a] 

[U] 

(31) 

tu] 

=  [a] 

[U] 

(32) 

With  the  assumed  velocity  distribution  within  the  element,  the 
"consistent  mass  matrix",  [M] ,  of  the  element  can  be  obtained 
easily  (25).  This  matrix  as  defined  by  Eq.  (12)  is 


(33) 


where,  the  various  symbols  appearing  in  the  above  matrix  have 
been  defined  previously  in  presenting  the  stiffness  matrices 
of  the  element. 

The  element  stiffness  and  mass  matrices  as  given  above 
(with  respect  to  the  element  local  coordinate  system)  can  now 
be  transformed  to  the  global  coordinate  system  by  utilizing 
the  transformation  matrix  represented  by  Eqs.  (22)  and  (23). 
The  matrix  operations  involved  in  such  a  transformation  are 
given  by  Eqs.  (18)  and  (19) . 
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IV.  THE  RECTANGULAR  PLATE  ELEMENT 


In  this  chapter  the  general  strain  energy  expression  for  a 
thin  rectangular  plate  subject  to  in-plane  and  bending  effects 
is  first  derived.  In  the  derivation,  the  interaction  between 
the  two  effects  is  also  taken  into  account.  The  stiffness  and 
mass  matrices  for  a  rectangular  plate  under  in-plane  action  and 
later  under  bending  action  are  derived  independently  by  assuming 
some  suitable  displacement  functions.  Following  that,  the  deri¬ 
vation  of  the  geometric  stiffness  matrix  is  given.  The  complete 
form  of  the  stiffness  and  mass  matrices  with  respect  to  both 
local  and  global  coordinate  systems  are  then  given. 

4.1  The  General  Strain  Energy  Expression 

In  order  to  derive  the  stiffness  properties  of  a  rectangular 
plate  in  which  coupling  between  in-plane  forces  and  bending  is 
considered  due  to  large  deformations,  it  is  convenient  to  obtain 
first  the  strain  energy  expression. 

The  middle  plane  of  a  plate  before  deformation  is  to  be  con¬ 
sidered  as  the  x-z  plane,  and  y  is  to  represent  the  axis  perpen¬ 
dicular  to  this  plane.  The  notation  used  in  subsequent  dis¬ 
cussions  is  similar  to  the  one  used  by  Timoshenko  (8). 

The  initial  stresses  and  strains  in  a  rectangular  element  of 

,  .  ,  ,  0  0  0,0  0  0 
the  plate  are  denoted  by  a  ,  a  ,  x  and  e  ,  e  ,  y  , 

X  Z  XZ  X  z  xz 

respectively.  It  is  assumed  that  stresses,  a  a  ®  x  ^ ,  have 

x  z  xz 

no  influence  on  subsequent  displacements  of  the  element  within 
the  initial  plane.  They  do,  however,  have  a  significant  effect 
on  displacements  of  the  element  out  of  its  initial  plane. 
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Starting  from  the  initial  position,  subsequent  deformation 
of  the  element  is  assumed  to  take  place  such  that 


0  a 

e  =  e  +  e 

X  XX 


0  ,  a 

e  =  e  +  e 
z  z  z 


(34) 


0  a 

Y  =  y  +  Y 
'xz  'xz  'xz 


^  a  a 

in  which  e  ,  e  ,  y  are  given  by  the  following  equations 
X  z  xz 

(8)  : 
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i  9u  2 
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3z 

2  l3z  ; 

-  y  ( — 2^) 

3z^ 

(35) 

a2u 


3u  3  u  3u  3u  . 

Yxza  =  +  +  (?x  }  "  2y{^¥)* 


In  the  above  equations  u  ,  u  and  u  represent  the  displacements 

x  y  z 

of  the  middle  surface  in  the  x,  y  and  z  directions,  respectively. 

The  strain  energy  expression  for  a  thin  rectangular  plate 
may  be  found  in  various  references,  such  as  Timoshenko  (8) ,  and 
is  in  the  form  of: 


V  =  i  ///  (as  +  a  e  +x  y  )dxdzdy  (36) 

2  xx  z  z  xz  'xz  2 

If  the  generalized  Hooke's  Law  for  an  isotropic  material  is  used, 
then 
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Here,  E  is  Young's  modulus;  v,  Poisson's  ratio  and  X  =  (l-v)/2. 
Substitution  of  Eqs.  (34),  (37)  into  Eq.  (36),  leads  to  the 
following  equation 


v  =  v0  +  vi  +  y2 


(38) 


where 


vo  ■  77r~27  S!S  "O  +  <S°>2  +  2v(ex0>  O 


2 (1-v  ) 


+  a(Yxz°)2]  dxdzdy 


(39) 


v  =  - ^-5-  Iff  (  (e  a) 2  +  (e_a)2  +  2v(e  a)(e  a) 

JL  o  /  n  . ,  ^  \  X  X  £* 


2 (1-V  ) 


+  X  (Yxj,3)2!  dxdzdy 


(40) 


V2  = 


(1-v") 


27  [ex°  C  +  V(exa  ez°  +  ex°  O 


+  ez°  eza  +  Ayxz°  Yxz3]  dxdzdy  (41) 


Note  that  Vq  is  a  constant,  and  since  its  second  order, 
partial  derivatives  are  zero  (see  Eq.  (11) ,  it  does  not  contribute 
to  the  stiffness  matrix.  It  can  therefore,  be  ignored  in  subsequent 
discussions . 

Suppose  the  following  assumptions  are  made: 

(1)  ux  =  ux(x,z) 

Uz  =  uz(x'z)  (42) 

Uy  =  Uy<X'Z) 
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(2)  The  cubic  and  higher  order  displacement  terms  are  ignored. 

The  first  assumption  is  a  general  assumption  for  solving 
the  plate  problem.  This  assumption  allows  the  elimination  of 
terms  containing  odd  powers  of  y  during  the  integration  process. 

The  second  assumption  limits  the  stiffness  matrix  to  that  of  a 
constant  matrix  (2) .  Substitution  of  Eqs.  (35)  into  Eq.  (40)  and 
utilizing  these  two  assumptions,  leads  to  the  following: 


V1  =  V1,P  +  V1,B 


(43) 


where 


V 


3u  2  3u  2 

Y  “  7  ^ 


1/P 


Eh  uuv  z  z  9u  9u  2 

27  "  [<3TT>  +  <inr>  +  x * +  <53r>> 


2 (1-V  ) 


9u 

+  2v(33r}  <3ir)]dxdz 


(44) 


V 


1/B 


32u  32u  2  92u  92uv 

iff  [(-/+-^)  -  2  (1-v)  ( — 2^  — 

z  3xZ  3zz  3xz  3zz 

32 

uv  2 

-  ( >  > 1 dxdz 


(45) 


In  Eqs.  (44)  and  (45)  the  constant,  h,  is  the  thickness  of 
the  plate  and  D  =  (Eh^) / (12  (1-v2) )  is  the  flexural  rigidity  of 
the  plate.  It  should  be  noted  that  the  two  terms  represented  by 
Eqs.  (44)  and  (45)  correspond  to  the  strain  energy  due  to  in-plane 
forces  and  bending  effects  respectively. 

Similarly,  substitution  of  Eqs.  (35)  into  Eq.  (41)  and  using 
the  same  assumptions  as  in  the  derivation  of  Eq.  (41) ,  the  follow¬ 
ing  expression  is  obtained. 
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(46) 


V2  =  V2  G  +  V2  1  +  V2  2 


where 


i  n  9u  •>  n  3u  2  n  3 u  9u 

4  w  to  U(^>  +  a,  (T^)  +  2t  °(^  +  ~-^)]dxdzdy  (47) 


n  9u  n  9u  n  9u  9u 

///  [ax°(9Fi)  +  +  Txz  <9lT  +  93r)^xdzdy  (48) 

92u  9  2u  9  2u 

-  Iff  y[ — ^  +  — 2^  +  2  9"x9 z^  dxdzdy  (49) 

3  x  9  z 


Note  that  V2  contains  only  terms  with  first  order  displace¬ 
ment  functions  and  hence  does  not  contribute  to  the  stiffness 
matrix.  The  integrand  of  V ^  2  an  oc^  function  of  y  and  hence 
the  value  of  the  integral  is  zero. 

Combining  results  from  Eqs.  (38)  to  (49)  the  following 
equation  is  obtained 

9u  2  3u  o  9u  9u  2 

v  =  - rT7//[<53r)  +  (*r>  +X(tt  +  kE) 

2  (1-v  ) 


3u  3u 

+  2v<93r>  W^dxdz 

n  9  2u  9  2u  2 


~  2  -2 

9  u  3  u 


+  5  "  + 

3x 

92u  2 

"  )]dxdZ 


- £)  -  2  (1-v)  ((— 5^  — 5- 

9zz  3x  9zz 


.  n  9u  2  n  9u  2  n  9u  9u 

+  2  //  [ax  (9^>  +  az  (9iZ)  +  2txz  (93EZ)  (nZ)]dxdz* 
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This  is  the  strain  energy  expression  due  to  the  coupling  of 
in-plane  and  bending  actions,  including  the  effect  of  large 
deformations.  The  first  term  on  the  right-hand  side  of  the  above 
equation  is  the  contribution  of  in-plane  action;  the  second  term 
represents  the  contribution  of  bending  action;  and  the  last  term 
is  the  contribution  due  to  large  deformations. 

Timoshenko  (8)  has  indicated  that  the  expression  represented 
by  Eq.  (50)  may  be  obtained  by  considering  the  strain  energy 
expression  separately  for  each  of  the  effects  due  to  in-plane, 
bending,  and  large  deformations.  These  results  are  then  added 
together  to  form  the  complete  strain  energy  expression.  The  authors 
believe  that  the  present  derivation  gives  a  better  insight  into  the 
problem. 

4 . 2  Coordinate  Systems  and  Degrees-of-Freedom 

A  Cartesian  coordinate  system  is  used  in  the  following 
derivation.  As  before,  let  the  global  coordinate  system  be  de¬ 
fined  as  x,y,z  (see  Fig.  2) .  In  general,  a  rectangular  plate 
element  may  be  oriented  in  any  arbitrary  direction  in  space, 
therefore,  coordinate  transformation  relations  have  to  be  estab¬ 
lished  as  explained  previously.  However,  for  simplicity  in  this 
discussion  let  its  middle  plane  be  considered  to  be  parallel  to 
the  X-Z  plane.  Let  the  origin  of  the  local  coordinate  system 
x,y,z  be  placed  at  one  of  the  corners  of  the  plate  such  that  the 
x-z  plane  coincides  with  the  middle  plane  of  the  plate.  In  a 
general  case,  at  each  corner  of  the  plate,  there  may  be  six 
possible  degrees-of-freedom.  For  convenience,  odd  numbers  are 
used  to  designate  those  pertaining  to  in-plane  action,  and  even 
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numbers  to  designate  those  pertaining  to  bending  action.  Fig.  2 
depicts  a  typical  plate  element  and  the  associated  element 
coordinates  corresponding  to  both  in-plane  and  bending  actions. 
It  may  be  observed  that,  in  a  general  case,  the  plate  element 
may  possess  a  total  of  24  degrees-of-f reedoms  corresponding  to 
both  in-plane  and  bending  actions.  Figs.  3  and  4  illustrate  the 
element  coordinates  corresponding  to  membrane  (in-plane)  and 
bending  actions,  respectively. 

In  what  follows,  stiffness  matrices  for  the  plate  elements 
are  derived  by  considering  in-plane  and  bending  actions  inde¬ 
pendently  of  each  other  on  the  basis  of  the  expressions  given  in 
the  previous  section.  The  combination  of  these  two  matrices 
together  with  the  so-called  geometrical  stiffness  matrix  gives 
the  complete  stiffness  matrix  of  the  structure.  However,  the 
complete  mass  matrix  contains  matrices  contributed  by  in-plane 
action  and  bending  action  only. 

4.3  Stiffness  Matrix  for  In-Plane  Action 

This  section  considers  the  stiffness  matrix  contributed  by 
in-plane  action.  Fig.  5  shows  the  plate  element  displacements 
corresponding  to  membrane  action.  The  corner  node  points  are 
labelled  in  a  clockwise  manner  as  i,j,k,i,,  so  that  node  i  coin¬ 
cides  with  the  origin  of  the  element  coordinate  system.  The 
x-dimension  of  the  rectangular  plate  (distance  from  node  i  to 
node  j ,  or  £  to  k)  is  designated  as  "a"  and  the  z-dimension 
(distance  from  node  j  to  node  k,  or  i  to  £)  is  designated  as  "b". 

Let  the  nondimensional  variables  £  and  r\  be  introduced,  so 

that  £  =  - 

8. 
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A  displacement  vector  U  is  defined  containing  the  membrane  dis¬ 
placements  at  the  corner  nodes  of  the  element,  so  that  (referring 
to  Fig.  5) 

U  —  (UjyU-^Uc^U^/Ug,  .  .  .  ,  t^^) 

The  displacements  in  the  x  and  z  directions  of  any  interior 

point  of  the  plate  are  referred  to,  as  before,  by  u  and  u  , 

respectively.  One  of  the  critical  steps  in  the  analysis  is  to 

find  suitable  representations  for  u  and  u  in  terms  of  the 

X  z 

nodal  displacement.  In  this  study,  the  following  displacement 
shape  functions  given  by  Weaver  et  al  (15)  are  used,  on  the 
basis  that  they  represent  shape  functions  that  are  compatible 
with  those  previously  employed  for  beam  elements. 


(l  -  O  ( 2n 3  -  3n2  +  l) 
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U1 

0 

u. 

3 

b  (l  -  £)  (n3  -  2n2  +  n) 

°5 
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U7 

0 

U9 

b?  (n3  -  2n2  +  n) 

U11 

-S  ( 2ti 3  -  3n2) 

U13 

0 

U15 

b£  (n3  -  n2) 

U17 

-U  -  ?)  (2n3  -  3n2) 

U19 

0 

U21 

b  (i  -  5)  (n3  -  n2) 

33_ 

31 


(51) 
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U9 

-a  (1  -  n)  (s3  -  d) 

U11 
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U13 

-n  (2£3  -  3C2) 

U15 
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U17 

0 

U,  „ 

19 

n  (2£3  -  3£2  +  1) 

U21 

-an  (d  -  2?2  +  £) 

_U23_ 

(52) 


These  two  displacement  functions  have  the  following  prop¬ 
erties  : 

(1)  For  an  arbitrary  set  of  nodal  displacements,  the  edges 
of  a  series  of  elements  describe  continuous  curves. 

(2)  For  an  arbitrary  set  of  nodal  displacements,  the  shear 
strain  at  each  corner  of  an  element  is  zero. 

(3)  They  are  compatible  with  displacement  functions  of 
beam  elements  described  by  Akkoush  et  al  (25) . 

Although  the  stiffness  matrix  [K]  of  the  plate  element  may 
be  obtained  by  partial  differentiation  of  the  strain  energy  ex¬ 
pression  given  by  Eg.  (50) ,  it  is  more  convenient  (15)  to  eval¬ 
uate  its  components  by  using  the  following  equation,  based  on  the 
principle  of  virtual  work: 
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(53) 


K..  =  h  //  [a.]T  [e.]  dxdz 

13  13 

where  the  index  j  represents  the  corresponding  strain  due  to 
a  unit  displacement  in  the  jth  direction  and  [a^]  represents 
the  stress  component  corresponding  to  the  strain  component  [e^] . 

Appendix  A  contains  the  various  strain  and  stress  components 
corresponding  to  the  displacement  functions  represented  by  Eqs. 
(51),  and  (52).  Thus,  the  stiffness  matrix  [ ]  contributed 
by  the  in-plane  action  with  displacement  functions  given  by 
Eqs.  (51)  and  (52)  can  be  obtained  by  the  substitution  into 
Eq.  (51)  of  the  stress  and  strain  components  given  in  Appendix  A. 
The  resulting  matrix  is  shown  below: 
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where,  the  various  symbols  appearing  in  the  above  equation  are 
defined  in  Table  1,  given  below  (a  and  b,  as  before,  pertain 
to  the  dimensions  of  the  rectangular  plate  element  and 
X  =  (l-v)/2,  where  v  is  Poisson's  ratio): 


TABLE  1  -  Definition  of  Symbols  Appearing  in  Equation  (54) 
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4.4  Mass  Matrix  for  In-Plane  Action 

If  the  same  displacement  functions  are  used  as  represented 
by  Eqs.  (51)  and  (52)  for  a  rectangular  plate  under  in-plane 
action,  then  the  corresponding  mass  matrix  [M^.]  can  be  obtained 
based  on  Eq.  (12) .  However,  for  the  plate  element  under  con¬ 
sideration,  the  matrix  [a]  in  Eq.  (12)  is  a  2  x  12  matrix  in 
which  the  first  and  second  rows  are  represented  by  the  first 
matrices  on  the  right-hand  side  of  Eq.  (51)  and  Eq.  (52) , 
respectively.  Since  the  procedure  involved  is  straight  forward 
(although  laborious) ,  only  the  final  result  is  presented  below: 
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(55) 
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where,  a  =  13/2520,  3  =  11/1260,  y  =  1/315,  6  =  1/1260,  p  =  3/140, 
and  X  =  1/840. 

4.5  Stiffness  Matrix  for  Bending  Action 

Many  displacement  functions  for  a  rectangular  plate  in  bending 
have  been  proposed.  Among  these,  the  one  which  satisfies  the 
deflection  and  slope  compatibility  has  been  chosen  in  this  study 
because  it  satisfies  not  only  compatibility  between  adjacent  plate 
elements  but  also  between  beam  and  plate  elements.  Such  a  dis¬ 
placement  function  representing  the  transverse  displacement  u^ 
of  the  plate  has  been  derived  and  verified  by  Bogner  et  al  (28) 
and  Przemieniecki  (27,29).  Fig.  6  illustrates  the  element  dis¬ 
placements  corresponding  to  bending  action. 
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For  a  plate  in  bending,  the  strain  components  have  the 
following  form: 
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xz 


=  -y 


32u 


3x 
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3  u 


9z' 


3  2u 


3x3z 


(57) 


The  matrix  [b]  given  in  Eq.  (2)  may  be  obtained  by  sub¬ 
stitution  of  Eq.  (56)  into  Eq.  (57) .  In  turn,  substitution  of 
matrix  [b]  into  Eq.  (9)  allows  evaluation  of  the  stiffness 
matrix  [Kg] ,  corresponding  to  bending  action. 

The  matrix  [Kg]  in  partitioned  form  is  represented  as 


[Kg] 


[K 

[K 


Bl- 


B2 ' 


Symmetric 

[kb33 


(58) 


The  submatrices  [Kg.^] ,  [Kg2] ,  and  [Kg3]  are  presented 
below  as  Eqs .  (60)  -  (62).  In  all  three  equations,  the  scalar 
multiplier  k  is  defined  as 


12 (1-v  ) ab 


(59) 


The  other  symbols  appearing  in  the  submatrices  are  defined  and 
compiled  in  Table  2. 
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TABLE  2  -  Definition  of  Symbols  Appearing  in  Equations  (60) -(62) 


i  =  1 

i  =  2 

i  =  3 

i  =  4 

i  =  5 

a  =  (b/a) 2 

w  =  (a/b)2 

6  =  l+5v 

e  =  l+60v 
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a .  = 
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78/35 

156/35 

xi  - 

1/50 

2/25 

6/25 

8/25 
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4.6  Mass  Matrix  for  Bending  Action 

In  generating  the  mass  matrix  for  the  plate  element  in  bending, 
if  the  same  displacement  function  given  by  Eq.  (56)  for  a  rectangular 
plate  in  bending  is  used,  then  the  corresponding  mass  matrix  [M  ] 

JD 

can  be  obtained  by  using  Eq.  (12).  The  matrix  [a]  in  Eq.  (12)  is 
a  1  x  12  matrix  represented  by  the  first  matrix  on  the  right-hand 
side  of  Eq.  (56)  .  The  resulting  matrix  [M_,]  is  shown  below. 
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where  the  various  symbols  appearing  in  Eq.  (63)  are  defined  as 
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and,  as  before,  a  and  b  refer  to  the  edge  dimensions  and  h  to  the 
thickness  of  the  plate,  and  p  represents  the  mass  density  of  the 


material 


4.7  The  Complete  Form  of  Stiffness  and  Mass  Matrices 

The  results  obtained  in  Sections  4.3  and  4.5  may  now  be 
incorporated  for  the  complete  form  of  the  local  stiffness 
matrix  of  a  rectangular  plate  with  coupling  between  in-plane 
and  bending  effects.  Let  [K^ ]  and  [Kg]  be  the  stiffness 
matrices  corresponding  to  in-plane  and  bending  actions,  re¬ 
spectively,  then,  the  force-displacement  relationship  for  the 
combined  effects  has  the  following  form: 

[S]  =  [K]  [U]  (64) 


where 

[S]  = 

and  [Sj]  ,  [Sg]  are  the  nodal  forces  corresponding  to  in-plane 
and  bending  actions,  respectively.  They  are  represented  by  the 
following  equations : 


[Sj] 


[Sg] 


(65) 


[Sj]  [s,  s3  s5  s?  s9  sn  s13  s15  s17  s19  s21  s23] 


T 


(66) 


[Sg]  -  [S2  S4  S6  Sg  S1(J  S12  S14  Slg  S18  S2Q  S22  S24]  .  (67) 


Similarly, 


[U] 


(68) 
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and  [UT] ,  [Uj  are  the  nodal  displacements  corresponding  to  in- 
X  D 

plane  and  bending  actions,  respectively.  They  have  the  following 
forms : 


IV  -  U1  u3  U5  U7  U9  U11  U13  U15  U17  U19  U21  U23! 


tuB]  =  u2  u4  u6  u8  u10  u12  u14  u16  u18  u20  u22  u24]  .  (70) 

Finally,  the  matrix  [K]  may  be  represented  as  follows: 


[K]  = 


[Kj]  [  0] 


[  0]  [KB] 


The  linear  orthogonal  transformation  [T]  is  introduced  for 
the  purpose  of  rearranging  the  forces  and  displacements  accord¬ 
ing  to  the  numbering  scheme  used  in  Fig.  2,  so  that: 


[T]  =  [  fT1]  [T2]  ] 

where  [T^] ,  [T2] ,  in  partitioned  form,  are  diagonal  matrices 

defined  as: 


[T1]  =  diag  ( [X1]  ,  [X1]  ,  [\±]  ,  [X]_]) 
[T2]  =  diag  ([A2J , [X2] , [X2] , [X2]) 


[A1]  = 


10  0 
0  0  0 
0  10 
0  0  0 
0  0  1 
0  0  0 
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(76) 


0  0  0 
10  0 
0  0  0 
0  10 
0  0  0 
0  0  1 

It  is  clear  from  the  definition  of  [T]  that 
[T] T  [T]  =  [T][T]T  =  [I] 

in  which  [T]  represents  the  transpose  of  [T] ,  and  [I] ,  the 
identity  matrix. 

Therefore,  if 

[S]  =  [K]  [U]  (77) 

then 

[T]  [S]  =  [T]  [K]  [U]  =  [T]  [K]  (  [T]  T  [T]  )  (U] 

=  (  [T]  [K]  [T]  T)  [T]  [U]  . 

Setting 

[S']  =  [T][S],  (78) 

[U>]  =  [T][U],  (79) 

and 

[K  *  ]  =  [T][K][T]T.  (80) 

Equation  (77)  is  thus  transformed  to 

[S']  =  [K '  ]  [U ' ]  (81) 


43 


which  is  the  complete  force-displacement  relationship  for  a 
rectangular  plate  with  coupling  between  in-plane  and  bending 
effects  in  accordance  with  the  numbering  scheme  indicated  in 
Fig.  2. 

The  complete  mass  matrix  for  a  rectangular  plate  with 
coupling  between  in-plane  and  bending  effects  again  according 
to  the  numbering  scheme  used  in  Fig.  2  may  be  obtained  in  a 
similar  way.  It  has  the  form  of: 

[M' ]  =  [T]  [M]  [T]T  ,  (82) 


where 

[M]  = 

The  matrices  [MT]  and  [FL,]  are  the  mass  matrices  corre- 

X  -D 

sponding  to  in-plane  and  bending  actions,  as  given  by  Eqs.  (55) 
and  (63) ,  respectively. 

4.8  Geometrical  Stiffness  Matrix 

It  has  been  indicated  above  that  the  stiffness  property  of 
a  rectangular  plate  is  contributed  by  three  terms  that  appear 
in  the  strain  energy  expression  for  the  plate.  Thus  far,  first 
and  second  term  contributions  have  been  dealt  with  in  the  deri¬ 
vations.  The  remaining  term,  i.e. ,  the  third  term,  is  the  one 
which  leads  to  the  geometrical  stiffness  matrix.  Since  this 


[MI]  [  0] 


[  0]  [Mfi] 


(83) 


term  involves  only  one  displacement  function  u 


Y ' 


it  is  necessary 
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to  assume  only  this  function.  The  derivation  of  this  geometrical 
stiffness  matrix  is  very  tedious,  consequently,  many  investigators 
have  tended  to  use  other  simpler  displacement  functions,  (e.g., 
see  Gallagher  et  al.  (30))  coupled  with  finer  partitioning  to  com¬ 
pensate  for  this  effect.  For  the  sake  of  being  consistent,  the 
same  displacement  function  as  shown  in  Eq.  (56)  is  used  here. 

The  strain  energy  contributed  by  large  deformations  which 
corresponds  to  the  third  term  on  the  right  of  Eq.  (50)  may  be 
rewritten  as  follows  (see  also  Eq.  (47)) 


3u 

(_Z) 

'  3  z  ' 


3u  3u 


)  ]  dxdz . 


The  derivation  of  the  geometrical  stiffness  matrix  may  be 
carried  out  in  three  parts,  each  of  which  corresponds  to  a  term 
in  Eq.  (84) .  Accordingly  then 


“V  -  V 


‘kg  1  +  V 

X 


tKGJ  + 


where  [KG  ] ,  [KG  ] ,  and  [KG  ]  are  the  parts  contributed  by  the 
x  z  xz 

first,  second,  and  third  terms,  respectively,  on  the  right  of 

Eq.  (84).  Utilizing  Eqs.  (56),  (11)  and  the  chain  rules  of 

differentiation,  the  matrices  [KG  ] ,  [KQ  ] ,  and  [KG  ]  may  be 

x  bz  xz 

obtained.  They  are  presented  below  as  Eqs.  (86)  -  (88) . 
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The  complete  geometrical  stiffness  matrix  may  now  be  obtained 

by  substitution  of  the  values  [K^  ] ,  [KQ  ] ,  and  [K^  ]  into  Eq.  (85)  . 

X  z  xz 

This  is  the  stiffness  property  contributed  by  the  large  deformation. 
Note  that  0  °,  0  and  t  0  are  treated  as  constants  during  the 
derivation  since  these  quantities  are  assumed  to  be  known  during 
the  current  state  of  deformation  of  the  structure. 

Since  the  geometrical  stiffness  matrix  has  contributions 
only  due  to  bending  action,  it  can  be  rearranged  easily  to 
correspond  to  the  numbering  scheme  indicated  in  Fig.  2. 

4.9  Transformation  to  the  Global _Coordinate  System 

In  order  to  assemble  individual  elements  to  form  a  complete 
structure,  the  local  properties  must  be  represented  with  respect 
to  the  global  coordinate  system.  It  is  seen  from  Eq.  (18)  that 
the  [X]  transformation  matrix  between  the  local  and  global  co¬ 
ordinates  has  to  be  determined  first.  The  stiffness  matrix  with 

respect  to  the  global  coordinate  system  may  then  be  obtained  by 

T 

evaluating  the  quantity  [X]  [K] [X] .  For  the  present  case,  let 

i,j,k,&  represent  the  four  corners  of  the  plate,  then 

[X]  =  diag  (  [ T i ]  ,  [Ti]  ,  (T  .]  ,  (T ^ ]  ,  [Tk]  ,  [TR]  ,  [T£]  ,  [T£]  )  (89) 


49 


where 


is  the  transformation  matrix  between  the  global  coordinate  system 
and  the  local  coordinate  system  at  vertex  "i".  This  transformation 
[T^]  is  actually  determined  by  the  direction  cosines  of  the  two 
edges  ij  and  iJl  (see  Figs.  5,6).  In  Eq.  (90),  x^,y^,z^  denote  the 
respective  coordinates  of  vertex  "i"  from  the  origin  and  ^ 
represent  the  distance  between  vertices  "i"  and  "j"  and  vertices 


"i"  and  respectively.  Moreover,  the  matrices  [T^],  [T^ ] ,  [T^] 

have  the  same  form  as  [T^]  except  that  they  are  determined  by  the 
direction  cosines  of  the  edges  ij  and  jk,  jk  and  ilk  and  i£  and  Jik, 
respectively,  rather  than  ij  and  i£. 

The  transformations  of  mass  and  geometrical  stiffness  matrices 
may  be  achieved  by  a  similar  procedure. 
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V.  PROCEDURE  OF  NUMERICAL  ANALYSIS  AND  THE  SOFTWARE  PACKAGE 


In  this  chapter,  first  the  general  procedure  is  outlined 
for  the  plastic  dynamic  analysis  of  a  structure  undergoing  large 
deformations.  This  is  followed  by  some  general  information  per¬ 
taining  to  the  software  package  that  has  been  developed  for  the 
analysis  of  the  class  of  structures  under  consideration. 

5.1  General  Procedure  of  Numerical  Analysis 

The  general  procedure  for  the  numerical  analysis  of  the 
problem  is  based  on  an  incremental  approach  which  allows  the 
computation  of  the  response  of  the  structure  at  discrete  time 
instances.  Within  each  interval  of  time  the  structure  is 
assumed  to  deform  as  a  linear  system  in  order  to  obtain  initial 
estimates  of  the  incremental  deformations.  Nonlinear  effects 
due  to  changes  in  the  geometry  of  the  structure  and  plastic 
deformations  are  then  taken  into  account  before  considering  a 
new  increment  of  time. 

It  is  assumed  that  the  response  of  the  structure  at  time 
t  =  t^  has  already  been  determined.  Considering  a  sufficiently 
small  increment  of  time,  At,  the  various  steps  for  computing 
the  response  of  the  structure  at  t2  =  t^  +  At  are  outlined  below: 

1.  Using  the  known  accelerations  at  t^,  the  displace¬ 
ments  at  t2  can  be  determined  based  on  a  numerical 
integration  algorithm  such  as  the  following 

Eq3  +.  =  tq]  +  At[q].  +  0.5(At)2  [q]  ,  (91) 

.  z2  rl  Z1  Z1 
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in  which  [q] ,  [q]  and  [q]  represent  the  generalized 

displacement,  velocity  and  acceleration  vectors, 
respectively.  The  displacement  increment  vector 
[Aq]  is  then  defined  as  follows 

[Aq]  =  [q]  .  -  tq]  .  (92) 

2  rl 

2.  The  incremental  element  displacement  vector,  [AU^] , 
with  respect  to  the  global  (generalized)  coordinate 
system,  is  determined  from  [Aq]  by  considering  the 
compatibility  between  generalized  and  element  dis¬ 
placements.  Then,  the  incremental  displacement 
vector,  [AU] ,  with  respect  to  the  local  (element) 
coordinate  system,  is  found  based  on  the  transformation 
relations  between  the  global  and  local  coordinate 
systems,  so  that 

[AU]  =  [X]  [ AU  ]  (93) 

where  [A]  refers  to  the  transformation  matrix  deter- 
rl 

mined  on  the  basis  of  the  deformed  configuration  of 
the  structure  at  t^. 

3.  Using  the  incremental  element  displacement  components, 
the  corresponding  incremental  forces  can  be  found 
easily  by  employing  the  force-deformation  properties 
of  the  elements.  However,  it  is  not  necessary  at  this 
time  to  compute  all  components  of  the  incremental  force 
vectors.  Thus,  in  the  case  of  a  beam  element,  only  the 


52 


force  increment  in  the  axial  direction  is  needed  in 


order  to  compute  the  geometric  stiffness  matrix  of 
the  element.  In  the  case  of  a  plate  element  (see 
Fig.  2)  the  components  pertaining  to 

the  forces  in  the  x-direction  and  S3  ,Sg  ,S21  in 
the  z-direction,  S4  'Sio,S16,S22  Pertain;*-n9  to  the 
moment  around  the  x-axis  and  Sg  S12'S18'S24  arounc^ 
the  z-axis  are  needed  to  compute  the  scalar  multi¬ 
plicative  factors  operating  on  the  elements  of  the 
geometric  stiffness  matrices.  Corresponding  approxi¬ 
mate  stress  components  a  0 ,  a  0 ,  t  ®  may  be  obtained 
by  dividing  the  forces  S  °,  S  °  and  S  0  by  the  thick- 
ness  of  the  plate.  These  forces  are  computed  based 
on  an  averaging  process  of  the  absolute  values  of  the 
force  components,  as  follows: 


,  1  ,  ' S7+S13 '  , 

+  2l  (  ^  + 
z 

,  1  # I S21+S15 I  , 

-  Tl~  (  2  + 

x 


sl+s19 


s3+s9 


1  ,  (S4+S10+S16+S22)  ,  (S'6+S12+S18+S24)  % 

L  L  (  4  4  ' 

X  z 


in  which  L  ,  L  are  the  dimensions  of  the  plate  in 
the  x  and  z  directions,  respectively.  The  sign  on 
the  right  hand  side  of  the  above  equations  is  taken 
as  positive  for  extension,  and  negative  for  compression. 
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Obviously,  this  procedure  of  obtaining  the  stress 

0  0  ,  0  .  . 
components  o  ,  a  ,  and  a  is  only  an  approximation; 

X  Z  X  2 

nevertheless,  it  has  been  utilized  by  other  investi¬ 
gators  (for  example,  see  Ref.  (30))  and  has  been  indi¬ 
cated  that  no  serious  error  will  occur  if  sufficiently 
large  numbers  of  idealized  elements  are  used  to  assemble 
the  structure. 

In  the  case  of  a  beam  element  the  geometric  stiff¬ 
ness  matrix  is  computed  based  on  Eq.  (30) .  The  cumula¬ 
tive  geometric  stiffness  matrix  (see  Eq.  (85)),  for  a 
plate  element  is  found  by  adding  together  the  three 
matrices  computed  from  Eqs .  (86)  -  (88).  The  geometric 

stiffness  matrix,  [KG] t  ,  computed  in  this  manner  for 
a  beam  or  plate  element  reflects  the  current  deformed 
configuration  of  the  structure  at  time  . 

4.  The  total  element  stiffness  matrix,  '  at 

is  then  formed  by  adding  the  elastic  and  geometric 
stiffness  matrices  for  the  element  under  consideration. 
Thus , 


[Ktot]  t. 


[K] 


+ 


[KG]t. 


(97) 


It  should  be  noted  that  although  the  elastic  stiff¬ 
ness  matrix  [K]  is  computed  based  on  the  element 
coordinate  system  and  the  small  deformation  (first 
order)  theory,  the  effect  of  the  large  deformation 
of  the  structure  enters  into  the  analysis  through 
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the  transformation  relations  pertaining  to  the  deformed 
and  undeformed  configurations.  Thus,  the  total  element 
stiffness  matrix  is  transformed  to  the  global  coordinate 
system  through  a  matrix  operation  represented  by  Eq.  (18) , 


i.e. 


r 


(Vt. 


[iWt. 


(98) 


where,  [Kt  ]  represents  the  total  element  stiffness 

matrix  with  respect  to  the  global  coordinate  system 

and  [A]  is  the  current  transformation  matrix  corre- 
t2 

sponding  to  the  deformed  configuration  of  the  structure 
at  t£. 

Now,  considering  the  entire  structure,  the  gener¬ 
alized  stiffness  matrix  for  the  entire  system 

is  obtained  in  the  usual  manner,  through  the  applica¬ 
tion  of  the  stiffness  method  (27) .  The  basic  process 
involved  is  that  of  considering  the  equilibrium  at  each 
nodal  point  of  the  forces  transmitted  by  all  elements 
which  are  incident  to  that  particular  node. 

5.  Although  some  parts  of  the  structure  may  undergo  plastic 
deformations,  in  order  to  obtain  estimates  of  the  ele¬ 
ment  forces,  initially  the  deformations  are  assumed  to 
be  entirely  elastic.  On  this  basis,  the  increments 
of  the  stress  resultants  acting  at  the  nodal  points  of 
the  element  can  be  found  to  be 

[ASI  -  tKtot]t2  1AU)  (99) 
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The  total  internal  forces  at  time  t  =  t2  are  then 
obtained  from 

ES]  =  [S],  +  [AS]  (100) 

rl 

6.  The  internal  stresses  computed  above  may  be  used  in 
connection  with  the  Mises-Hencky  yield  criterion  (22) 
to  determine  whether  the  element  under  consideration 
undergoes  elastic  or  inelastic  deformations.  In  terms 
of  the  stress  resultants  this  criterion  is  expressed  as 

$  (Si,S2, _ ,S y...)  =  Y  (101) 

where  S^,S2,...,  represent  the  normalized  form  of 
the  element  forces  S^,S2,...,  and  Y  denotes  the 
"yield  value"  which  may  change  through  straining  (for 
simplicity,  the  subscript  t2  is  omitted  here  and  in 
subsequent  discussions,  except  where  necessary). 

7.  If  the  yield  function  $  is  less  than  Y  the  element 
under  consideration  is  deforming  elastically.  However, 
if  the  yield  function  is  greater  than  Y  or  equal  to  Y 
and  its  rate  of  change  is  positive,  then  plastic  defor¬ 
mations  are  taking  place.  If  the  material  of  the 
element  is  assumed  to  have  isotropic  strain  hardening 
properties,  subsequent  yielding  of  the  element  will 
occur  whenever  $  attains  a  value  equal  to  or  greater 
than  the  "current"  value  of  Y.  However,  the  consider¬ 
ation  of  kinematic  hardening  properties  is  physically 
more  realistic  and  accounts  for  the  Bauschinger  effect. 
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If  all  elements  are  found  to  behave  elastically 
within  the  current  time  increment,  then  the  acceler¬ 
ations  can  be  found  using  Eq.  (20) ,  but  omitting  the 
term  [F°] .  Then  returning  to  Step  1,  the  above  process 
can  be  repeated  for  the  next  increment  of  time. 


8.  The  stress  resultants  used  in  forming  $  have  been  found 
using  stress  increments  based  on  elastic  considerations. 
In  the  case  of  plastic  deformations,  these  stresses, 
being  only  crude  approximations  to  the  actual  values, 
do  not  satisfy  the  plastic  flow  relations.  Therefore, 
for  a  plastically  deforming  material  the  correct  values 
of  the  stress  resultants  have  to  be  approximated  more 
closely  by  using  the  procedure  explained  below. 

As  discussed  previously,  it  is  assumed  that  any 
typical  displacement  increment  AU^  may  be  expressed  as 
the  sum  of  an  elastic  and  a  plastic  part,  i.e., 


AUj  =  AU j  +  AU?  (102) 

The  increments  of  the  elastic  components  [AUe]  of 
each  element  are  related  to  the  internal  force 
increments  through  the  inverse  of  the  total  element 
stiffness  matrix 

[AUe]  =  [Ktot]_1  [AS]  (103) 

fc2 
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9.  The  j-th  plastic  component  of  the  displacement  incre¬ 
ment  vector  may  be  found  using  the  generalized  form  of 
the  Reuss-Mises  plastic  flow  rule,  so  that 

AU?  =  A  —  (104) 

3  as . 

J 

where  X  is  a  constant  of  proportionality  which  usually 
changes  through  straining. 

This  approach  is  applicable  to  the  general  case 
of  a  stress-strain  curve  as  shown  in  Fig.  7.  However, 
in  order  to  simplify  the  numerical  computations  a  piece- 
wise  linear  stress  strain  curve  may  be  used.  For  ex¬ 
ample,  the  trilinear  curve  indicated  in  Fig.  8  is  a  fre¬ 
quently  used  approximation.  Let  H'  represent  the  inelastic 
loading  curve  and  S  the  slope  of  the  unloading  curve.  The 
initial  estimates  of  the  force  components  may  be  computed  as 

(AS  .)  -  H'  AU.  (105) 

3  est  3 

10.  The  elastic  and  plastic  components  of  the  displacements 
are  then  found  to  be 

e  (ASi ) pqf 

AuJ  =  - 5~—  ■■  (106) 

J  ^ 

AU?  =  AU  .  -  AUe  (107) 

HDD 

11.  Using  the  current  element  stiffness  matrix  (see 

Eq,  (97)),  the  final  values  of  the  force  components 
(stress  resultants)  may  be  found  from 
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(108) 


[AS]f  =  [Ktot][AU6] 

12.  The  above  procedure  is  repeated  for  all  other  elements 
deforming  plastically.  The  vector  [AU°]  =  [AUP]  is 
then  formed  and  used  to  find  the  increments  of  the 
equivalent  generalized  nodal  forces  due  to  plastic 
deformations ,  i . e . , 

[ AF° ]  =  [ktot][q°]  =  [X] T (Ktot] [AU°]  (109) 

Also, 

[F° ]  .  =  [F° ]  +  [ AF° ]  (110) 

*1 

13.  The  acceleration  vector  [q]  at  t  =  can  now  be  found 
based  on  the  governing  equation  of  motion  for  the 
entire  system  as  given  by  Eq.  (20) .  Also  the  velocity 

vector  [q] .  is  determined  by  numerical  integration. 

2 

14 .  The  entire  procedure  explained  above  may  be  repeated 
for  a  new  time  increment  At,  and  in  fact  for  any 
prescribed  number  of  time  increments.  Thus,  the  history 
of  internal  forces  and  displacements  of  the  structure 
for  a  given  time  interval  can  be  determined. 

5.2  General  Remarks 

It  should  be  emphasized  that  the  computation  of  the  plastic 
displacement  components  as  indicated  in  Steps  9-11  of  the  above 
procedure  represents  a  somewhat  simple  and  perhaps  crude  method 
of  calculation  of  these  components .  The  main  intent  has  been  to 
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minimize  to  the  extent  possible  the  size  of  a  fairly  complex 
problem  in  order  to  realize  solutions  that  are  economically 
feasible  from  the  computational  viewpoint.  However,  for  com¬ 
parison  purposes  a  more  sophisticated  method  of  computation  of 
plastic  deformations  has  also  been  utilized.  This  alternative 
method  is  based  on  an  iterative  scheme  aimed  at  examining  the 
accuracy  of  estimated  quantities  at  successive  steps  of  iteration. 

In  the  iterative  scheme,  initial  estimates  of  the  incre¬ 
mental  element  forces  are  utilized  in  order  to  compute  the  plastic 
and  elastic  components  of  the  displacement  increments  from 
Eqs.  (104)  or  (106)  -  (107) .  An  initial  value  of  the  effective 
strain  is  then  calculated  by  utilizing  well  known  relations, 
as  given  for  example,  by  Hill  (22)  .  New  estimates  of  the  incre¬ 
mental  element  forces  are  now  found  from  Eq.  (108)  and  are  used 
to  calculate  a  second  set  of  plastic  and  elastic  displacement 
components.  A  new  value  of  the  effective  strain  is  then  computed 
and  compared  with  the  previously  determined  value.  If  the  dif¬ 
ference  between  the  two  values  is  within  a  specified  tolerance, 
then  the  iterative  process  is  terminated.  Otherwise,  the  new 
estimates  of  the  incremental  element  forces  are  used  towards  a 
new  iterative  step. 

Comparison  of  the  two  methods  has  been  made  in  the  case  of 
a  simple  structure  which  was  analyzed  by  means  of  the  computer 
program  discussed  in  Section  5.3  of  this  report.  Results  indicate 
that  there  is  no  significant  difference  in  the  two  methods.  How¬ 
ever,  before  a  final  conclusion  can  be  reached  more  extensive 
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numerical  experimentation  should  be  undertaken. 

It  should  also  be  pointed  out  that  the  same  general  procedure 
outlined  in  Section  5.1  can  be  used  in  the  case  of  relatively 
small  deformations  in  the  plastic  range  by  neglecting  the  effect 
of  the  geometric  stiffness  matrix,  i.e. ,  not  taking  into  account 
the  matrix  [KQ] .  This  is  equivalent  to  formulating  the  equilibrium 
equations  with  respect  to  the  original  (undeformed)  configuration 
of  the  structure.  In  this  manner,  the  need  of  recomputing  the 
total  stiffness,  mass,  and  transformation  matrices  at  specified 
time  steps  is  eliminated. 

Simplified  forms  of  the  general  procedure  in  Section  5.1  may 
be  used  to  treat  the  case  of  the  elastic  deformations  of  the 
structure.  Thus,  the  procedures  involved  in  the  elastic  static 
and  elastic  dynamic  analyses  of  a  structure  may  be  deduced  easily 
from  the  general  procedure  by  following  through  those  steps  that 
are  appropriate  to  each  particular  type  of  analysis.  It  should 
be  emphasized,  however,  that  in  the  case  of  the  static  large 
deformation  analysis,  increments  in  the  loading  system  rather 
than  time  have  to  be  considered.  For  each  increment  in  the  load 
vector  the  corresponding  increment  in  the  displacement  vector  is 
computed  and  added  to  the  previously  accumulated  displacements. 

The  current  total  stiffness  and  transformation  matrices  are  then 
determined  prior  to  considering  a  new  load  increment.  This 
process  is  continued  until  the  loads  attain  their  final  values. 
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5.3  Plastic  Behavior  of  In-Plane  Stressed  Plate  Element 

The  implementation  of  the  general  procedure  in  Section  5.1 
to  the  case  of  a  plate  element  undergoing  plastic  deformations 
poses  some  practical  difficulties.  If  the  most  general  type  of 
plate  element  with  in-plane,  bending  and  twisting  action  is  used, 
then  the  element  possesses  24  degrees-of-f reedom  (see  Fig.  2), 
as  indicated  previously.  Thus,  if  a  plate  subassembly  of  an 
actual  structure  is  modelled  by  several  plate  elements  of  this 
type,  the  number  of  degrees-of-f reedom  for  the  structure  under 
consideration  would  increase  rapidly.  This  may  not  be  as  sig¬ 
nificant  in  the  case  of  an  elastic  analysis  of  the  structure, 
but  tends  to  be  critical  when  attempting  a  plastic  dynamic 
analysis  of  the  system,  due  to  excessive  computational  times 
involved.  In  such  a  case,  a  method  of  reducing  the  size  of  the 
problem  is  highly  desirable.  In  this  study,  a  technique  has 
been  developed  for  treating  a  rectangular  plate  element  with 
in-plane  action  as  a  single  degree-of-f reedom  system  only  when 
plastic  deformations  take  place. 

Consider  the  plate  element  shown  in  Fig.  9,  with  8  degrees- 
of-freedom  at  the  nodal  points  corresponding  to  in-plane  action 
with  no  rotational  displacements  at  the  corners.  The  degrees- 
of-freedom  are  labelled  according  to  the  numbering  scheme  indicated 
in  Fig.  2,  and  the  corresponding  stress  resultants  (generalized 
stresses)  are  designated  as  ,  S3,  S?,  Sg  ,  S13,  S15 ,  Sig  and  S21. 
Considering  the  stress  resultants  acting  on  each  edge  of  the  plate 
an  average  shear  stress  can  be  found  as  an  approximation  to  the 
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stresses  acting  over  the  area  of  the  edge.  Thus,  referring  to 
Fig.  9 
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An  average  shear  stress  for  each  pair  of  parallel  edges  is 
then  found  using  the  stresses  computed  from  Eqs.  (Ill)  and  utiliz¬ 
ing  the  usual  sign  convention  in  the  theory  of  elasticity.  Denot¬ 
ing  by  tab  and  tcd  the  average  shear  stresses  acting  on  the  edges 
parallel  to  the  x  and  z  axes,  respectively 


„  _^ItaI  +  Itb 

AB  +  2 


(112) 


TCD  - 


Tc ' + ' TD I 


(113) 


where,  the  minus  sign  in  Eq.  (112)  is  used  when  ta  is  positive, 
while  the  plus  sign  in  Eq.  (113)  corresponds  to  the  case  of  tc 
being  positive. 

Finally,  an  average  shear  stress  for  the  entire  plate 
element  is  computed  as  follows 


TAV  =  ± 


tab ' + ' tcd ' 


(114) 
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The  sign  in  Eq.  (114)  is  chosen  to  concide  with  the  sign  of 
either  xAfi  or  tCD-  This  is  based  on  whether  I TAB ! > I Tcp L  or 
I  x  ,  I  <  I  x  I.  In  the  former  case,  the  sign  of  x„TT  is  taken  to 
be  the  same  as  that  of  and  in  the  latter  case  it  is  chosen 

to  be  the  same  as  that  of  x^D.  In  implementing  the  above  pro¬ 
cedure  in  a  computational  algorithm  several  intermediate  checks 
are  performed  to  verify  that  the  plate  is  in  equilibrium  in 
both  the  x  and  z  directions  within  specified  tolerances. 

The  average  shear  stress  found  from  Eq.  (114)  for  each 

plate  element  is  then  used  in  connection  with  the  procedure  in 

Section  5.1  to  determine  whether  plastic  deformations  are  taking 

place  in  the  element.  Also,  in  the  computer  program,  is 

automatically  checked  against  the  ultimate  stress  of  the  material 

(see  Fig.  8)  to  determine  whether  the  element  has  failed  due  to 

excessive  straining.  If  the  ultimate  stress  is  not  exceeded 

but  plastic  deformations  have  been  found  to  take  place,  then  an 

initial  estimate  of  the  actual  shear  stress,  x  . ,  in  the  element 

est 

is  computed  from 


T 


est 


S 


Y 

-) 


H' 


(115) 


where.  x ,  refers  to  the  known  shear  stress  in  the  element  at  the 


1 

end  of  time  t^,  and  as  before,  H' 
the  inelastic  loading  and  elastic 
and  Y  is  the  current  yield  value, 
term  on  the  right  hand  side  of  Eq 
the  sign  of  x,  . 

t  -i 


and  S  pertain  to  the  slope  of 
unloading  curve,  respectively, 
The  sign  in  front  of  the  second 
(115)  is  chosen  to  agree  with 
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The  equivalent  element  nodal  force  increment  vector,  [AS0] , 

due  to  plastic  deformations  can  now  be  found  based  on  the  usual 

procedure  in  the  application  of  the  initial  strain  method  (20) , 

i.e.,  referring  to  the  jth  equivalent  element  force  increment 

AS  ? 

3 

AS?  =  (s.)  (lA-y-T-e-St-)  (116) 

3  3  t2  AV 

where,  the  (S.)  's  are  again  calculated  from  Eqs.  (99)  and  (100). 

3  t2 

The  transformation  of  the  equivalent  element  force  increment 
vectors  of  all  plate  elements  to  the  global  coordinate  system 
and  consideration  of  the  equilibrium  conditions  at  each  nodal 
point  yields  the  vector  [AF°] .  The  usual  procedure  outlined  in 
Section  5.1  is  then  followed.  As  explained  in  the  previous 
section,  an  iterative  scheme  can  also  be  used  in  order  to  improve 
on  the  values  of  the  estimated  stresses. 

5.4  The  Software  Package 

A  general  purpose  computer  program  has  been  developed  for  the 
analysis  of  frame  and  plate  structures.  This  program  represents 
an  extension  of  a  previously  reported  (31)  computer  program  known 
as  the  GWU-FAP  (George  Washington  University  -  Frame  Analysis 
Program)  which  deals  with  the  analysis  of  rigid  frame  structures. 
Extensions  to  GWU-FAP  have  been  accomplished  in  order  to  incor¬ 
porate  a  rectangular  thin  plate  element  and  to  develop  the  capa¬ 
bility  for  second  order  (large  deformation)  analysis  that  accounts 
for  geometric  nonlinearities.  In  addition,  the  plastic  dynamic 
analysis  branch  of  the  program  has  been  augmented  in  order  to 
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allow  for  the  analysis  of  structures  consisting  of  in-plane 
stressed  rectangular  plate  elements  (see  Section  5.3).  A  simpli¬ 
fied  form  of  the  flow  chart  for  the  modified  GWU-FAP  is  presented 
in  Appendix  B. 

The  modified  GWU-FAP  is  a  general  purpose  program  which 
enables  the  user  to  perform  the  following  types  of  analysis  for 
both  frame  and/or  plate  structures  (with  certain  exceptions  in 
the  case  of  a  plastic  dynamic  analysis  as  explained  below) : 

(a)  Elastic  static  small  deformation  (first  order)  analysis 

(b)  Elastic  static  large  deformation  (second  order)  analysis 

which  considers  the  effect  of  geometric  nonlinearities 

(c)  Free  vibration  analysis  leading  to  the  natural  frequencies 

and  modes  of  vibration,  using  the  consistent  mass 
matrix  of  the  structure 

(d)  Elastic  dynamic  first  order  analysis 

(e)  Elastic  dynamic  second  order  analysis 

(f)  Plastic  dynamic  first  order  analysis  of  structures 

consisting  of  beam  or  in-plane  stressed  rec¬ 
tangular  plate  elements. 

When  performing  any  of  the  above  types  of  analysis  for  a 
rigid  (skeletal)  frame,  the  structure  is  modelled  as  an  assembly 
of  plane  or  three-dimensional  beam  elements,  depending  on  the 
loading  and  deformation  states  of  the  structure.  When  con¬ 
sidering  a  structure  which  contains  rectangular  plate  elements, 
different  types  of  plate  action  may  be  distinguished.  In  the 
present  study,  five  possible  modes  of  plate  action  have  been 
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considered  and  are  classified  according  to  the  following  scheme 
Class  1  -  In-plane  action  without  rotational  displacement 
at  the  nodes. 

This  type  of  element  possesses  8  degrees-of- 
freedom  corresponding  to  element  coordinate 
directions  1,3,7,9,13,15,19  and  21  indicated 
in  Fig.  2  (see  also  Fig.  9) . 

Class  2  -  In-plane  action  with  rotational  displacements 
at  the  nodes. 

This  type  of  element  possesses  12  degrees-of- 
freedom  corresponding  to  element  coordinate 
directions  1,3,5,7,9,11,13,15,17,19,21  and  23 
indicated  in  Fig.  2  (see  also  Fig.  3) . 

Class  3  -  Bending  action  only. 

This  type  of  plate  element  also  has  12  degrees- 
of-freedom  represented  by  coordinate  directions 

2.4.6.8.10.12.14.16.18.20.22  and  24  in  Fig.  2 
(see  also  Fig.  4) . 

Class  4  -  Bending  and  in-plane  action  without  rotational 
displacements  at  the  nodes. 

This  type  of  plate  element  possesses  20  degrees- 
of-freedom  represented  by  element  coordinate 
directions  1,2,3,4,6,7,8,9,10,12,13,14,15,16,18, 

19.20.21.22  and  24  in  Fig.  2. 
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Class  5  -  Bending  and  in-plane  action  with  rotational 


displacements  at  the  nodes. 

This  type  of  plate  element  is  the  most  general 
type  used  in  the  present  study  and  possesses 
24  degrees-of-f reedom  as  indicated  in  Fig.  2. 

Any  of  the  above  classes  of  rectangular  plate  elements 
may  be  used  when  performing  an  analysis  of  type  (a)  through 
(e)  for  a  structure  consisting  of  plate  elements  or  a  combina¬ 
tion  of  beam  and  plate  elements.  On  the  other  hand,  in  the 
case  of  a  structure  containing  rectangular  plate  elements ,  a 
plastic  dynamic  analysis  corresponding  to  type  (f)  can  be  per¬ 
formed  only  with  Class  1  plate  elements. 

The  computer  program  has  been  designed  so  as  to  minimize 
the  amount  of  input  needed  for  any  of  the  above  types  of 
analysis.  Most  of  the  input  information  to  be  supplied  by  the 
user  of  the  program  pertains  to  the  type  of  the  structure,  its 
geometry,  its  material,  cross-sectional  and  inertial  properties, 
the  type  and  duration  of  loading  and  the  type  of  analysis  de¬ 
sired.  The  input  format  for  some  typical  structures  corre¬ 
sponding  to  different  types  of  analysis  is  given  in  Appendix  C. 

The  computer  program  is  written  in  Fortran  IV  language  for 
use  on  an  IBM-OS/360  digital  computer.  It  has  been  made  oper¬ 
ational  on  both  the  APL  Model  91  and  the  George  Washington  Uni¬ 
versity  Model  50  Computer  System. 
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VI.  COMPUTER  ANALYSIS  OF  TYPICAL  STRUCTURES-SUMMARY  OF 

NUMERICAL  RESULTS 


The  software  package  discussed  in  the  previous  chapter  is 
used  here  to  obtain  solutions  to  typical  problems.  For  this 
purpose,  different  types  of  structures  are  analyzed  and  the 
results  of  the  numerical  solutions  are  summarized  in  this  chapter. 
Several  of  the  solutions  presented  in  the  initial  part  of  the 
chapter  pertain  to  simple  structures  and  serve  the  purpose  of 
checking  the  various  branches  of  the  computer  program  by  comparing 
with  known  solutions. 

6.1  First  Order  Static  Analysis  of  a  Plate 

One  of  the  simplest  type  of  structures  considered  in  this 
study  for  the  purpose  of  testing  the  computer  program  consists 
of  a  simply  supported  plate  with  a  concentrated  load  at  the  center. 
Both  rectangular  and  square  plates  with  dimensions  (length  x  width 
x  thickness)  40  x  20  x  0.50  in  and  20  x  20  x  0.50  in,  respectively, 
are  used.  The  other  constants  are: 

Young's  Modulus:  E  =  30  x  10^  lb/in^ 

Poisson's  Ratio:  v  =  0.3 
Concentrated  load:  P  =  10,000  lb. 

The  type  of  plate  element  used  corresponds  to  Class  3, 
described  in  Section  5.4.  As  mentioned  previously,  this  type 
of  plate  element  is  subjected  only  to  bending  action,  and 
possesses  3  degrees-of-f reedom  at  each  node  (corner)  point,  or, 
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a  total  of  12  degrees-of-freedom.  In  order  to  compare  the  accuracy 
of  the  numerical  results,  solutions  are  obtained  for  the  rectangular 
plate  by  subdividing  it  into  4,  16,  and  32  plate  elements;  similarly, 
solutions  are  obtained  for  the  square  plate  by  subdividing  it  into 
4  and  16  elements. 

Table  3  shows  the  comparison  of  the  center  deflections  of  the 
plate  obtained  by  the  finite  element  method  and  the  exact  method 
as  given  in  Timoshenko  (8) : 


TABLE  3  -  STATIC  ANALYSIS  OF  SIMPLY  SUPPORTED  PLATE— COMPARISON 
OF  CENTRAL  DEFLECTION  BETWEEN  FINITE  ELEMENT 
AND  EXACT  METHODS 


Displacements  (in) 

Square  Plate 

Rectangular  Plate 

a  =  20  in 

a  =  20  in 

b  =  40  in 

4  Elements  .11767 

.16920 

16  Elements  .12564 

.17706 

32  Elements 

.18332 

Timoshenko  (6)  .1351 

.1923 

(exact) 

It  may  be  observed  that  in  both  cases  the  solutions  appear 
to  converge  to  the  exact  values  as  the  number  of  elements  is 
increased. 

6.2  Free  Vibration  Analysis  of  a  Simply  Supported  Plate 

A  free  vibration  solution  is  obtained  for  a  simply  supported 
square  plate  with  the  following  physical  input: 
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Young's  Modulus:  E  =  30  x  10^  lb/in2 

Dimensions  of  plate:  a  x  b  x  h  =  20  x  20  x  0 . 50  in 

Poisson's  Ratio:  v  =  0.3 

2 

Density  of  plate:  p  =  0.001  lb-sec  /in. 

Again,  a  Class  3  plate  element  is  used  in  modelling  the 
structure  which  is  subdivided  into  4  and  16  elements. 

The  theoretical  results  used  for  comparison  are  obtained 
from  Volterra  (32)  .  The  eigenvalues,  wmn*  and  the  corre¬ 
sponding  eigenvectors,  Wmn(x,y),  for  the  problem  under  con¬ 
sideration  are  given  as  follows: 


w 


ran 


(117) 


W 

mn 


(x,y) 


sin  (525*)  sin  (2£*> 


(118) 


where 


8  =  Eh2/ ( 1 2 p (1-v2)  )  , 


is  a  numerical  constant. 

Table  4  shows  the  comparison  of  the  fundamental  frequencies 
between  results  obtained  by  the  use  of  Eq.  (117)  and  the  finite 
element  method  as  applied  in  this  study. 


TABLE  4  -  FREE  VIBRATION  OF  SIMPLY  SUPPORTED  SQUARE  PLATE- - 
COMPARISON  OF  FUNDAMENTAL  FREQUENCIES  BETWEEN 
FINITE  ELEMENT  AND  EXACT  METHODS 


4  elements 

16  elements 

Volterra (22) - (exact) 

Fundamental  Frequency 
Corresponding  Period 

.  217630E0cyc/sec 

.459496E-02  sec. 

.  212534E03cyc/sec 

.470513E-02  sec. 

. 205913E03cyc/sec 

. 485642E-02  sec. 
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It  may  be  seen  that  even  with  the  use  of  a  small  number  of  ele¬ 
ments  the  finite  element  solutions  are  fairly  close  to  the  exact 
solution. 

Table  5  shows  the  comparison  of  the  fundamental  normal  mode 
shape  between  the  finite  element  (16  elements)  and  the  exact 
method. 


TABLE  5  -  FREE  VIBRATION  OF  SIMPLY  SUPPORTED  SQUARE  PLATE— COMPARISON 
OF  FUNDAMENTAL  NORMAL  MODE  SHAPES  BETWEEN 
FINITE  ELEMENT  AND  EXACT  METHODS 


Volterra  (32)  (exact) 

16  Elements 

(39  degrees-of-freedom) 

-0 . 536411E-01 

-0 . 478539E-04 

-0 . 549793E-01 

-0 . 665452E-07 

-0. 758360E-01 

0.340921E  00 

-0 . 777525E-01 

0.340488E  00 

-0.535734E-01 

0. 338058E-04 

-0 . 549791E-01 

0.110401E-06 

0 . 536411E-01 

-0.536749E-01 

0 . 549792E-01 

-0 . 549792E-01 

0 . 241449E  00 

-0. 758359E-01 

0 . 240762E  00 

-0 . 777516E-01 

-0 . 379180E-01 

0. 535734E-01 

-0 . 388762E-01 

0 . 549787E-01 

0 . 379180E-01 

0.241144E  00 

0.388764E-01 

0 . 240760E  00 

0 . 341353E  00 

0. 379659E-01 

0 . 340489E  00 

0 . 388763E-01 

-0 . 536072E-01 

0. 378701E-01 

-0 . 549791E-01 

0 . 388761E-01 

-0. 338485E-04 

0.340921E  00 

-0. 503340E-07 

0 . 340487E  00 

0 . 241144E  00 

0.536749E-01 

0 . 240762E  00 

0 . 549793E-01 

-0. 378701E-01 

-0.338058E-04 

-0. 388761E-01 

-0 . 116132E-06 

-0. 379659E-01 

0.240840E  00 

-0 . 388764E-01 

0 . 240760E  00 

-0. 536410E-01 

0. 379179E-01 

-0 . 549790E-01 

0 . 388762E-01 

0 . 758360E-01 

-0 . 379179E-01 

0. 777519E-01 

-0 . 388762E-01 

0 . 341353E  00 

-0.535733E-01 

0.340488E  00 

-0 . 549784E-01 

0.338485E-04 

0.536410E-01 

0 . 144227E-07 

0 . 549785E-01 

0 . 536072E-01 

0.758359E-01 

0 . 549793E-01 

0 . 777512E-01 

0 . 482593E  00 

0 . 535733E-01 

0.481525E  00 

0 . 549783E-01 

0 . 478539E-04 

0 . 443147E-07 
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Comparison  between  the  two  methods  indicates  that  in  most  cases 
the  agreement  is  good,  except  when  the  amplitude  of  vibration  is 
very  small.  This  is  attributed  to  rounding  and  truncation  errors 
in  the  computational  process. 

A  comparison  of  higher  frequencies  corresponding  to  the  two 
methods  is  shown  in  Table  6 . 

TABLE  6  -  FREE  VIBRATION  OF  SIMPLY  SUPPORTED  SQUARE  PLATE --COMPARISON 
OF  HIGHER  FREQUENCIES  OF  VIBRATION  BETWEEN 
FINITE  ELEMENT  AND  EXACT  METHODS 


4  Elements 

16  Elements 

Volterra  (22) 

(m,n)  (exact) 

.217630E03/sec. 

. 212534E03/sec . 

mm 

.  205913E03/sec . 

. 614070E03/sec. 

. 530795E03/sec . 

HI 

. 514782E03/sec. 

. 614070E03/sec. 

. 530796E03/sec. 

.  514  782E03/sec. 

.  123725E04/sec . 

. 870523E03/sec . 

(2,2) 

.  82365lE03/sec. 

. 1067 59E0 4/sec . 

(1,3) 

.  102956E04/sec . 

.106759E04/sec. 

(3,1) 

.  102956E04/sec . 

.131904E04/sec. 

•142567E04/sec. 

(2,3) 

. 133843E04/sec. 

. 142567E04/sec. 

(3,2) 

. 133843E04/sec. 

It  is  worth  noting  that,  even  with  only  16  plate  elements,  the 
comparison  of  higher  frequencies  obtained  in  this  study  with 
the  frequencies  given  in  Volterra  (32)  is  in  general  favorable 
even  for  the  higher  frequencies  of  vibration. 

6.3  Second  Order  Static  Analysis  of  a  Square  Plate 

The  second  order  static  analysis  of  a  square  plate  under  a 
concentrated  load  at  the  center  is  studied  in  this  section. 
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Two  types  of  boundary  conditions  are  used  for  this  structure: 

(1)  all  four  edges  clamped  and  (2)  all  four  edges  simply  supported 

The  plate  dimensions  are  30  x  30  x  0.25  in,  Young's  modulus  and 

6  2 

Poisson's  ratio  are  taken  as  30  x  10  lb/in  and  0.3,  respectively 
For  the  purpose  of  comparison  with  the  work  of  other  investi¬ 
gators,  the  type  of  plate  element  chosen  is  of  Class  4  (see 
Section  5.4).  This  type  of  plate  element  is  subjected  to  bending 
and  in-plane  action  without  rotational  degrees  of  freedom  at  the 
nodes.  The  ultimate  load  is  50,000  lb.,  and  16  elements  are  used 
to  idealize  the  square  plate.  Numerical  results  corresponding 
to  the  clamped  plate  solution  are  compared  with  results  obtained 
by  Brebbia  and  Connor  (12) ,  and  Adotte  (33) .  Table  7  below, 
shows  the  similarities  and  differences  among  the  three  methods. 

TABLE  7  -  COMPARISON  OF  METHODS  USED  IN  SOLVING  THE  NONLINEAR  PROBLEM 
OF  A  SQUARE  PLATE  CLAMPED  AT  FOUR  EDGES  WITH 
A  CONCENTRATED  LOAD  AT  THE  CENTER 


Adotte 

Brebbia  and  Connor 

Present  Study 

Method 

Finite  difference 
and  experimental 

Finite  element 

Finite  element 

Number  of 
Elements 

— 

36 

16 

Degrees-of- 
freedom  at 
each  corner 

5 

5 

In-plane 

displacement 

functions 

_  „ , 

simpler 

more  complex 

Numerical 

procedure 

- - - 

used  linearized 
equations  for  a 
limited  number  of 
load  steps,  then 
applied  correc¬ 
tions  based  on  a 
Newton-Raphson 
_ iterative  method 

100  piecewise 
linear  in¬ 
cremental 
steps 
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Figure  10  shows  the  results  obtained  by  the  three  different 
methods.  It  is  seen  that  Adotte ' s  curve  is  enveloped  by  the 
curve  obtained  in  the  present  study  and  that  of  Brebbia  and 
Connor,  with  the  former  below  and  the  latter  above  it.  All 
three  results  are,  however,  very  close,  and  the  characteristic 
of  the  second  order  analysis  in  which  the  structure  becomes 
stiffer  as  the  load  increases  is  clearly  evident. 

Figure  11  shows  the  results  of  both  the  clamped  plate  and 
simply  supported  plate  plotted  on  a  different  scale  with  the 
coordinate  axes  switched  around.  Both  curves  show  the  same 
characteristic  of  increasing  stiffness  of  the  structure  with 
increasing  load. 

6.4  First  and  Second  Order  Static  and  Elastic  Dynamic 

Analysis  of  a  Beam-Plate  Assemblage 

In  this  section  is  presented  the  analysis  of  a  beam-plate 
assemblage  including  the  first  and  second  order  static,  the  free 
vibration,  and  the  first  and  second  order  elastic  dynamic 
analysis.  The  structure  consists  of  a  horizontal  square  plate 
elastically  supported  by  four  edge  beams  which,  in  turn,  are 
supported  in  the  vertical  direction  by  four  edge  beams.  The 
plate  is  subjected  to  a  uniformly  distributed  load. 

The  following  numerical  values  are  used  for  the  problem 
under  consideration: 

Plate  dimensions:  30  x  30  x  0.25  in 

Poisson's  ratio:  v  =  0.3 

1=0 

x 

El  =  aD,  where  a  is  the  length  of  square  plate  and 
D  the  plate  stiffness. 
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Both  a  total  of  4  and  16  Class  3  plate  elements  are  used 

to  idealize  the  square  plate  in  order  to  compare  results  with 

solutions  given  by  other  investigators,  such  as,  Timoshenko  (8) 

and  Zienkiewicz  (13).  For  this  purpose,  initially  a  uniformly 

2 

distributed  load  q  =  1  lb/in  is  used.  A  comparison  of  the 
results  is  given  in  Table  8. 

TABLE  8  -  FIRST  ORDER  STATIC  ANALYSIS  OF  A  SQUARE  PLATE  WITH  FOUR 
EDGE  BEAMS --COMPARISON  OF  RESULTS 
WITH  KNOWN  SOLUTIONS 


Center  Point 
Deflection 

Deflection  at  mid-point 
of  edge 

4  Elements 

0.1274 

0.0626 

16  Elements 

0.1584 

0.0713 

v  =  0.30 

Timoskenko  (exact) 

0.1639 

not  available 

v  =  0.25 

Zienkiewicz 

0.1639 

0.0697 

36  Elements 

v  =  0.30 

It  is  seen  from  the  above  table  that  the  16  element 
solution  obtained  in  this  study  agrees  fairly  closely  with 
the  solutions  given  by  Timoshenko  and  Zienkiewicz. 

.  2 

6.4.1  Second  Order  Static  Analysis  -  A  uniform  load  of  10  lb/m 
is  now  applied  to  the  plate  in  order  to  study  the  effect  of 
geometric  nonlinearities.  For  this  purpose  a  second  order  analysis 
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is  performed  using  a  total  of  100  incremental  steps.  Fig.  12 
shows  the  force-deformation  curve  for  the  center  of  the  plate 
and  for  the  mid-point  of  the  edges.  For  comparison  purposes, 
also  the  first  order  force-deformation  curve  is  plotted.  It 
may  be  observed  again  that  the  structure  tends  to  become  stiffer 
as  the  intensity  of  the  load  is  increased. 

6.4.2  Free  Vibration  Analysis  -  The  free  vibration  analysis  of 
the  structure  under  consideration  yields  the  following  results: 

Lowest  Frequency  =  0.766  cy c/sec. 

Largest  Period  =  1.31  sec. 

Highest  Frequency  =  410.6  cyc/sec. 

Smallest  Period  =  0.0024  sec. 

6.4.3  First  and  Second  Order  Elastic  Dynamic  Analysis  -  The 
frequencies  and  periods  from  the  free  vibration  analysis  may 
be  used  as  a  guideline  for  the  numerical  input  to  the  elastic 
dynamic  analysis.  Thus,  a  time  increment  of  0.001  seconds 

and  a  final  time  of  2.7  seconds  is  used  for  the  numerical  inte¬ 
gration  process  as  applied  to  the  first  order  analysis.  In 
the  case  of  the  second  order  dynamic  analysis,  a  smaller  time 
increment  has  to  be  used  since  the  structure  becomes  stiffer 
as  the  deformations  increase,  due  to  stretching  of  the  middle 
surface  of  the  plate.  It  has  been  found  that  a  time  increment 
of  0.0001  seconds  satisfies  the  stability  criterion  for  the 
numerical  integration  algorithm  and  yields  satisfactory  results. 

For  the  first  order  dynamic  analysis,  uniformly  distributed 
loads  of  intensity  1  and  10  lb/in  are  used.  In  the  case  of  the 
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.  2 
second  order  analysis  load  intensities  of  1,2,5  and  10  lb/in 

are  utilized  in  conjunction  with  an  incremental  procedure.  For 

comparison  purposes  the  time-history  curves  from  both  analyses 

are  plotted  on  the  same  figure  (Fig.  13) . 

Curves  and  S^q  in  Fig.  13  correspond  to  the  first  order 

2 

dynamic  analysis  for  load  intensities  of  1  and  10  lb/in  ,  re¬ 
spectively.  Within  a  time  interval  of  2.4  seconds,  two  nearly 
identical  complete  cycles  are  seen  for  each  curve,  each  possessing 
a  period  of  oscillations  of  about  1.2  seconds,  but,  of  course, 
differing  in  amplitude  by  a  factor  of  approximately  10.  The 
corresponding  static  displacements  can  be  found  from  Fig.  12  and 
although  not  plotted  on  Fig.  13,  they  represent  horizontal 
equilibrium  lines  for  the  dynamic  oscillations. 

Curves  and  L1Q  in  Fig.  13  represent  results  obtained 

based  on  a  second  order  analysis  for  similar  load  intensities 

2 

of  1  and  10  lb/in  ,  respectively.  Those  two  curves  clearly 

demonstrate  the  difference  between  the  first  and  second  order 

analyses.  Within  the  same  time  interval  of  1.2  seconds,  approx - 

2 

imately  four  complete  cycles  appear  for  the  1  Ib/in  loading  and 

2 

twelve  complete  cycles  for  the  10  lb/in  .  The  maximum  ampli¬ 
tudes  are  much  smaller  than  those  of  the  first  order  analysis. 

This  is  because  the  plate  becomes  stiffer  due  to  the  stretching 
of  the  middle  surface.  It  is  also  seen  that  initially,  the 
results  for  both  first  and  second  order  analyses  are  quite  close 
(as  may  be  observed  by  comparing  curve  S^q  with  L-^q  and  S.^  with 
L^) .  However,  after  some  time  has  elapsed  (for  a  load  intensity 
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2 

of  10  lb/in  ,  this  corresponds  to  about  0.2  seconds,  and  in  the 
2 

case  of  1  lb/in  ,  to  about  0.4  seconds),  the  solutions  of  the 

first  and  second  order  dynamic  analyses  become  markedly  different 

Figures  14  shows  the  comparison  of  second  order  dynamic 

analyses  for  uniformly  distributed  loads  of  different  intensities 

Curves  L^,  ,  L^,  and  L^g  represent  the  time  history  of  the 

center  deflection  of  the  square  plate  under  uniformly  distributed 

2 

loads  of  1,  2,  5,  and  10  lb/in  ,  respectively. 

These  curves  show  the  following  characteristics: 

(1)  All  curves  are  periodic,  or  at  least  tend  to  be 
periodic,  after  some  initial  time  has  elapsed,  since  no  damping 
has  been  introduced  into  the  analysis. 

(2)  The  period  of  oscillations  becomes  shorter  as  the 
load  increases.  There  are  about  five  complete  cycles  in  L1  but 
about  fifteen  or  more  complete  cycles  in  L^g  for  the  same  time 
interval.  Again,  this  reveals  the  fact  that  the  plate  is  much 
stiffer  and  vibrates  much  faster  as  the  load  increases. 

(3)  The  amplitudes  of  oscillations  increase  but  in  a 
nonlinear  manner  due  to  the  nonlinear  stiffness  properties  of 
the  structure. 

6.5  Square  Plate  Subjected  to  a  Lateral  Load  and 
In-Plane  Compressive  Forces 

In  this  section  is  discussed  both  the  elastic  static  and 
dynamic  analyses  of  a  simply  supported  square  plate  under  a 
lateral  concentrated  load  at  center  and  uniform  compressive 
forces  in  the  x-direction. 
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Up  to  now,  only  one  direction  of  loading  has  been  considered. 
When  a  plate  is  under  the  effect  of  both  a  lateral  load  and  a 
uniform  compressive  load  acting  in  its  middle  plane,  the  work  done 
by  both  the  transverse  load  and  the  compressive  forces  has  to 
be  considered  in  formulating  the  solution  to  the  problem.  This 
problem  may  be  solved  by  using  the  analytical  expressions  given  in 
Chapter  IV  since  the  work  done  by  the  compressive  forces  is  equi¬ 
valent  to  the  last  term  on  the  right  hand  side  of  Eq.  (50) .  In 
order  to  compare  results  with  Timoshenko  (16)  ,  the  case  of  a 
simply  supported  square  plate  under  a  concentrated  load  at  center 
and  uniformly  compressed  in  the  x-direction  is  considered.  The 
deflection  surface  of  a  simply  supported  rectangular  plate  as 
given  by  Timoshenko  is : 
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and,  P  represents  the  concentrated  force  at  the  point  (£  ft)  / 

Nx  is  the  compressive  force  in  the  x-direction. 

It  is  seen  that  when  N  =  0 ,  the  above  equation  reduces  to 

X 

the  solution  of  a  simply  supported  rectangular  plate  under  a 
lateral  concentrated  load.  In  order  to  compare  results,  the 
concentrated  load,  P,  is  kept  constant  at  1000  lb.,  acting  at 
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the  center  of  the  square  plate,  whereas  N  is  varied.  In  the 
present  solution  by  the  finite  element  method,  in  most  cases  a 
total  of  16  Class  4  (see  Section  5.4)  plate  elements  are  used 
to  idealize  the  square  plate.  The  numerical  constants  used  in 

the  first  order  static  analysis  are: 

6  2 

Young's  Modulus  E  =  30  x  10  lb/in  . 

Poisson's  ratio  v  =  0.3. 

Dimensions  of  plate  =  30  x  30  x  0.25  in. 

5 

Flexural  rigidity,  D  =  0.429258  x  10  lb-in. 

As  long  as  the  denominator  on  the  right  hand  side  of  Eq.  (120) 
does  not  approach  zero,  the  double  series  on  the  right  hand  side 
of  Eq.  (119)  converges  rapidly.  Hence  the  solution  obtained  by  the 
finite  element  method  is  compared  only  with  the  first  term  of  the 
series  using  the  coefficient  A^,  and  the  results  of  the  compari¬ 
son  are  shown  in  Table  9. 

TABLE  9  -  FIRST  ORDER  STATIC  ANALYSIS  OF  A  SIMPLY  SUPPORTED  SQUARE  PLATE 
UNDER  A  1000  LB.  LOAD  AT  CENTER  AND  UNIFORM  COMPRESSIVE 
FORCES  IN  THE  X-DIRECTION- -COMPARISON  OF  CENTRAL 
DEFLECTION  BETWEEN  FINITE  ELEMENT 
AND  EXACT  METHODS 


Deflection  -  Inches 

F.E. 

Method 

Timoshenko  (24) 
(exact) 

N  (lb/in) 

A 

4  Elements 

16  Elements 

0 

—  —  mm  — 

.22615 

.21489 

100. 

— 

.23694 

.22694 

200. 

— 

.24892 

.24041 

400. 

— 

.27727 

.27280 

600. 

— 

.31358 

.31528 

800. 

— 

.36180 

.37342 

1000. 

.39909 

.42901 

.45787 

1200. 

— 

.52926 

.59166 

1400. 

— 

.69523 

.83592 

1600. 

---- 

1.02366 

1.42366 

1800. 

---- 

1.98451 

4.7953 

1900. 

~  — - 

00 
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It  may  be  observed  that,  in  general,  the  comparison  is 
fairly  good  for  N  less  than  1200  lb/in.  However,  when  N  is 

X  X 

greater  than  1200  lb/in,  increases  very  rapidly  due  to  the 

fact  that  its  denominator  tends  to  become  very  small.  This 
is  because  the  critical  value  of  N  is  between  1800  to  1900  lb/in. 

The  same  structure  is  now  analyzed  by  performing  a  second 
order  static  analysis  to  determine  the  center  deflection  of  the 
plate.  The  loads  are  increased  to  their  final  values  by  apply¬ 
ing  10  equal  increments.  Table  10  shows  a  comparison  of  the 
accuracy  between  a  4  and  16  element  idealization  of  the  square 
plate,  with  P  =  1000  lb. ,  and  N  =  1000  lb/in. 

X 

TABLE  10  -  SECOND  ORDER  STATIC  ANALYSIS  OF  A  SIMPLY  SUPPORTED  SQUARE  PLATE 
WITH  P  =  1000  LB  AT  THE  CENTER  AND  UNIFORM  COMPRESSIVE 
FORCES  Nx  =  1000  LB/ IN— COMPARISON  OF 

CENTRAL  DEFLECTION  FOR  4  AND  16 
ELEMENT  IDEALIZATIONS 


Deflection  -  Inches 

Step 

16  Elements 

4  Elements 

1 

.02369 

.02221 

2 

.04807 

.04503 

3 

.07188 

.06729 

4 

.09397 

.08793 

5 

.11372 

.10642 

6 

.13107 

.12272 

7 

.14624 

.13704 

8 

.15957 

.14968 

9 

.17134 

.16089 

10 

.18182 

.17091 
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The  values  in  the  last  line  of  Table  10  represent  the  deflections 
at  the  center  of  the  plate  when  the  loads  have  attained  their  full 
values.  It  is  seen  that  the  difference  between  the  two  idealizations 
is  about  6%. 

The  effect  of  geometric  nonlinearities  can  best  be  studied 
by  varying  both  the  lateral  load  and  in-plane  compressive  forces 
applied  in  the  x-direction.  For  this  purpose,  values  of  P  =  500, 
1,000,  2,000,  10,000  lb.  and  N  =  0,  1,000,  2,000,  3,000  and 
4,000  lb/in.  are  used.  The  results  of  the  investigation  are 
summarized  in  Table  11. 

TABLE  11  -  SECOND  ORDER  STATIC  ANALYSIS  OF  A  SIMPLY  SUPPORTED  SQUARE  PLATE 
WITH  VARYING  LATERAL  AND  COMPRESSIVE  FORCES— 

COMPARISON  OF  CENTRAL  DEFLECTIONS 


Deflection  -  Inches 

N 

X 

P 

0 

1000 

2000 

3000 

4000 

500 

.0932 

.1118 

.1364 

.1644 

.1904 

1000 

.1542 

.1709 

.1891 

.2079 

.2264 

2000 

.2312 

.  2440 

.25695 

.2700 

.28292 

10000 

•48755 

.49481 

.50235 

_ 

.51026 

.51865 

It  may  be  observed  from  the  above  table  that  for  smaller 
lateral  concentrated  loads,  the  effect  of  varying  N  on  the 
deflection  of  the  plate  is  much  more  pronounced.  For  large 
values  of  P,  large  deformations  of  the  plate  are  accompanied  by 
stretching  of  the  middle  surface.  This,  in  turn,  nearly  counter¬ 
acts  the  compressive  effect  of  the  in-plane  forces. 


83 


An  elastic  dynamic  analysis  has  also  been  performed  for  the 
same  structure.  Initially,  a  free  vibration  solution  is  obtained 
which  indicates  that  the  largest  period  is  about  0.7  seconds. 
Therefore,  a  time  interval  of  1.4  seconds  is  used  for  the  numeri¬ 
cal  solution  process  in  the  forced  vibration  analysis.  The  time 
increment  used  is  0.001  seconds.  In  Fig.  15  are  shown  the  time- 
history  curves  corresponding  to  the  first  order  dynamic  analysis. 

Curves  D..  ,  S  are  the  dynamic  and  static  displacement  curves 
%  wx 

for  P  =  2000,  N  =  1000;  D  ,  S  are  the  dynamic  and  static  dis- 
X  o  o 

placement  curves  for  P  =  2000,  Nx  =  0.  It  is  seen  that  has 

X 

a  much  greater  amplitude  and  slightly  longer  period  of  oscillations 

than  Dq.  In  fact,  the  fundamental  period  obtained  from  the  free 

vibration  analysis  no  longer  matches  the  period  of  oscillations 

corresponding  to  response  curve  D  ,  due  to  the  existence  of 

Nx 

the  in-plane  compressive  forces. 

6.6  First  Order  Elastic  and  Plastic  Dynamic  Analysis  of 
a  Plate  Structure 

In  several  of  the  numerical  solutions  presented  in  this 
chapter,  it  has  been  assumed  that  the  stresses  and  strains  in 
every  part  of  the  structure  remain  elastic  even  when  large 
deformations  are  taking  place.  In  such  cases  it  is  more  realistic 
to  take  into  account  plasticity  effects  which  are  brought  about 
by  the  yielding  of  the  material.  To  accomplish  this  the  "Plastic 
Dynamic"  branch  of  the  computer  program  discussed  in  Section  5.4 
is  used  to  solve  some  typical  problems. 
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One  of  the  simple  structures  considered  for  this  purpose 
consists  of  a  square  plate  which  is  fixed  along  one  of  its  edges 
and  free  along  the  others  (see  Fig.  16) .  The  dimensions  of  the 
plate  are  30  x  30  x  0.50  in.  The  structure  is  discretized  using 
4  plate  elements  corresponding  to. Class  1,  as  discussed  in 
Section  5.4.  This  type  of  plate  element  can  sustain  only  in¬ 
plane  action  without  twist  at  the  corner  nodes.  As  shown  in 
Fig.  16,  three  concentrated  in-plane  loads  of  magnitude  F  =  100  kips 
are  applied  in  the  x-direction  at  the  nodes  C,  D  and  E  on  the 
free  edge  of  the  plate  parallel  to  the  fixed  support.  For  com¬ 
parison  purposes,  first  a  static  analysis  of  the  structure  is  per¬ 
formed.  The  results  indicate  that  the  maximum  displacement  due 
to  the  applied  loads  is  0.146  in.  and  occurs  in  the  x-direction 
at  the  nodes  C  and  E.  The  corresponding  displacement  at  node  D 
is  0.137  in.  Likewise,  in  order  to  obtain  an  estimate  of  the  time 
increment  of  integration  to  be  used  in  the  dynamic  analysis  of 
the  structure,  a  free  vibration  analysis  is  first  performed.  The 
results  of  this  analysis  are 

Lowest  Frequency  =  25.88  cyc/sec. 

Largest  Period  =  0.0386  sec. 

Highest  Frequency  =  233.7  cyc/sec 
Smallest  Period  =  0.0043  sec. 

Based  on  the  above  results,  an  elastic  and  plastic  dynamic 
analysis  of  the  structure  is  performed.  The  same  concentrated 
loads  used  in  the  static  analysis  are  now  applied  to  the  struc¬ 
ture  dynamically  and  are  assumed  to  remain  constant  with  time. 
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A  time  increment  of  0.0002  sec.  is  used  for  the  numerical  inte¬ 
gration  process.  In  the  case  of  the  plastic  dynamic  analysis  a 
trilinear  stress-strain  curve  is  used.  Referring  to  Fig.  8,  the 
values  of  the  various  parameters  defining  the  curve  are: 


Param. 

1: 

Stress  = 

20,000 

lb/in2 

Param. 

2: 

Slope  2  = 

11.50 

x  106  lb/in2 

Param. 

3: 

Stress  a2  = 

33,000 

lb/in2 

Param, 

4: 

Slope  3  = 

2.50  x 

106  lb/in2 

Param. 

5: 

Ultimate  Stress  = 

130,000  lb/in2 

Param. 

6: 

Yield  Point 

Stress 

=  33,000  lb/in 

It  should  be  pointed  out  that  initially  a  more  realistic  value 
2 

of  65,000  lb/in  was  used  for  the  ultimate  stress  of  the  material. 
However,  it  was  found  that  the  average  stresses  in  some  elements 
exceeded  the  ultimate  stress.  This,  in  turn,  led  to  an  automatic 
interruption  of  the  computer  analysis,  as  explained  in  Section  5.3. 

The  time-history  curve  of  the  x-displacement  at  node  C  for 
both  the  elastic  and  plastic  dynamic  analysis  is  plotted  in 
Fig.  17.  A  comparison  of  the  two  curves  indicates  that  in  the 
case  of  plastic  deformations  there  is  an  increase  of  about  5% 
in  the  maximum  displacement  and  a  slight  increase  in  the  period 
of  oscillations.  However,  the  main  difference  between  the  results 
of  the  two  analyses  is  that  plastic  deformations  are  accompanied 
by  an  noticeable  permanent  set,  as  may  be  seen  from  Fig.  17. 

6.7  Elastic  and  Plastic  Dynamic  Analysis  of  a  Plane  Frame 

The  structure  considered  in  this  case  consists  of  a  plane 
frame,  as  indicated  in  Fig.  18(a),  that  was  previously  analyzed 
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in  Ref.  (25) .  The  frame  is  made  up  of  standard  4  inch  diameter 

2 

steel  pipe  whose  cross-sectional  area  A  =  3.174  m  and  the 

4 

moment  of  inertia  I  =  7.233  in  .  Furthermore,  the  modulus  of 

6  2 

elasticity  E  is  taken  to  be  29  x  10  Ib/in  ,  and  a  value  of 
2 

46,000  lb/in  is  used  for  the  yield  value  of  the  material.  For 
simplicity  in  the  computation  of  plastic  deformations  a  bilinear 
stress-strain  curve  is  used,  with  the  inelastic  branch  having  a 
slope  of  1/10  of  that  of  the  elastic  branch.  The  loading  consists 
of  two  dynamically  applied  forces,  F  =  20  kips  and  F„  =  5  kips 
having  a  rise  time  of  0.20  seconds,  as  shown  in  Fig.  18(b). 

The  time-history  of  deflection  at  the  point  of  application 
of  the  vertical  load  is  shown  in  Fig.  19.  Also  included  in  the 
figure  is  the  corresponding  deflection  if  the  loads  are  applied 
statically  and  the  structure  is  assumed  to  respond  as  a  linear 
system.  The  three  time  history  curves  shown  in  the  figure  corre¬ 
spond  to  the  elastic  dynamic  response  with  or  without  geometric 
nonlinearity  effects,  and  the  plastic  dynamic  response  of  the 
structure.  As  may  be  observed,  the  dynamic  character  of  the  load 
results  in  a  dynamic  load  factor  of  approximately  2.0  for  the 
case  where  the  effects  of  geometric  nonlinearities  are  neglected 
(first  order  analysis) .  However,  when  such  effects  are  retained 
(second  order  analysis)  the  load  factor  increases  to  approximately 
2.4  (20%  increase).  In  the  case  of  inelastic  deformations,  the 
dynamic  load  factor  becomes  much  larger.  This  is  due  to  the 
large  magnitude  of  the  applied  loads  which  causes  initiation  of 
yielding  at  an  early  stage  of  the  response  (i.e.,  at  approximately 
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0.16  sec.)  and  thus  leads  to  considerable  inelastic  deformations. 


6.8  Elastic  and  Plastic  Dynamic  Analysis  of  a  Three-Dimensional 

Structure  Consisting  of  In-Plane  Stressed  Plate  Elements 

The  static,  elastic  and  plastic  dynamic  analysis  of 
the  structure  shown  in  Fig.  20  is  studied  in  this  section.  The 
structure  is  composed  of  aluminum  plates  of  uniform  thickness. 

It  consists  of  three  vertical  walls  (parallel  to  the  y  -  z  plane) 
at  6  ft.  apart  and  three  horizontal  floor  decks  at  2  ft.  apart. 

The  plate  thicknesses  chosen  are  1/8  in.  for  the  vertical  walls 
and  upper  deck  and  1/16  in.  for  the  lower  two  decks.  The  struc¬ 
ture  is  supported  rigidly  along  the  bottom  edges  of  the  vertical 
walls  and  is  restrained  from  movement  in  the  x-direction. 

A  discrete  model  of  the  structure  is  obtained  by  inserting 
a  node  at  each  corner  point  corresponding  to  the  interconnection 
between  vertical  and  horizontal  members.  In  this  manner,  each 
vertical  wall  is  subdivided  into  three  and  each  deck  into  two 
plate  elements.  These  plate  elements  are  of  the  type  that  can 
sustain  only  in-plane  stresses  (see  Sections  5.3  and  5.4)  with 
no  twisting  action.  Thus,  at  each  node  point  only  a  translational 
movement  in  the  y  and  z-directions  is  allowed,  and  therefore  the 
number  of  degrees  of  freedom  for  the  entire  structure  is  36. 

The  structure  is  subjected  to  a  concentrated  load  of 
F  =  900  kips  applied  at  node  16  and  acting  in  the  positive  z- 
direction.  Initially,  the  load  is  considered  to  act  statically 
in  order  to  determine  the  static  displacement  distribution  of  the 
structure.  Numerical  results  indicate  that,  as  expected,  the 
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maximum  displacement  occurs  in  the  direction  of  the  applied  load 
and  is  equal  to  0.  887  in.  Also,  a  free  vibration  analysis  of 
the  structure  has  been  performed  leading  to  the  following  results: 

Lowest  Frequency  =  236.2  cyc/sec. 

Largest  Period  =  0.00423  sec. 

Highest  Frequency  =  4540  cyc/sec. 

Smallest  Peroid  =  0.00022  sec. 

It  may  be  observed  that  the  three-dimensional  character  of  the 
structure  results  in  a  fairly  stiff  system  and  corresponding  high 
frequencies  of  vibration.  Although  not  plotted  here,  a  study 
of  the  fundamental  mode  of  vibration  indicates  that  the  structure 
deforms  in  a  manner  analogous  to  a  cantilever  beam.  Thus,  as 
far  as  the  z-direction  is  concerned,  all  nodal  points  move  to¬ 
gether  in  the  same  direction.  However,  the  movement  in  the 
y-direction  is  such  that  while  all  points  in  the  frontal  plane 
(parallel  to  the  x-y  plane)  move  in  one  direction,  all  points 
in  the  rear  plane  (x-y)  plane  move  in  the  opposite  direction. 

The  same  concentrated  load  used  in  the  static  analysis  is 
now  applied  to  the  structure  dynamically  and  is  assumed  to  remain 
constant  with  time.  The  elastic  dynamic  response  of  the  structure 
is  first  obtained  by  assuming  the  structure  to  act  as  a  linear 
elastic  system.  The  time-history  curve  of  the  displacement  in 
the  direction  of  the  applied  load  is  shown  in  Fig.  20.  Some 
secondary  oscillatory  effects  are  observed  in  the  plot  due  to  the 
interaction  between  inertial  effects  in  the  two  perpendicular 
directions  and  the  superposition  of  higher  modes  of  vibration. 
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For  the  plastic  dynamic  analysis  a  trilinear  stress-strain 
curve  is  used  as  represented  by  Fig.  8.  The  following  values  of 
the  parameters  are  used  in  this  case: 


Param. 

1: 

Stress  a ^  = 

6,530 

lb/in2 

Param. 

2: 

Slope  2  = 

3.87  x 

106  lb/in2 

Param. 

3: 

Stress  = 

20,800 

lb/in2 

Param. 

4: 

Slope  3  = 

0.03  x 

106lb/in2 

Param. 

5: 

Ultimate  Stress 

=  30,000  lb/in2 

Param. 

6: 

Yield  Point 

Stress 

=  20,800  lb /in2 

The  same  magnitude  of  the  load,  i.e.,  F  =  900  kips,  is 
used  in  the  plastic  dynamic  analysis  of  the  structure.  The 
time-history  curve  of  the  displacement  in  the  direction  of 
the  load  is  plotted  also  in  Figure  21.  The  same  oscillatory 
behavior  is  observed  as  the  case  of  the  elastic  dynamic 
analysis,  however,  the  displacements  are  larger.  Also, 
plastic  action  is  accompanied  by  relatively  large  permanent 
deformations,  indicating  possible  damage  to  the  structure. 


90 


VII.  SUMMARY  AND  CONCLUSIONS 


A  procedure  of  analysis  has  been  presented  for  determining 
the  plastic  dynamic  response  of  structural  systems  consisting  of 
beam  and  rectangular  plate  elements.  The  analysis  takes  into 
account  both  the  geometric  and  material  nonlinearities. 

The  general  approach  to  the  problem  is  based  on  the  finite 
element  method  using  displacement  interpolation  functions.  The 
strain  energy  expressions  for  both  the  beam  and  plate  elements 
are  obtained  and  are  used  to  generate  the  stiffness  and  mass 
matrices  of  the  elements.  Also,  the  geometric  stiffness  matrices 
are  derived  that  account  for  the  effect  of  geometric  nonlinearities. 
Plastic  deformations  are  taken  into  account  by  means  of  an  incre¬ 
mental  theory  of  plasticity  coupled  with  the  concept  of  initial 
strain.  The  governing  dynamical  equation  of  the  system  is 
written  based  on  Hamilton's  Principle. 

A  general  purpose  computer  program  has  been  constructed 
based  on  the  analytical  procedures  developed  in  this  study. 

Although  the  program  is  primarily  intended  for  the  solution  of 
plastic  dynamic  problems,  it  can  also  be  used  to  perform  an 
elastic  static  or  elastic  dynamic  analysis  with  or  without 
geometric  nonlinearities,  as  well  as  a  free  vibration  analysis. 

The  program  is  used  in  this  study  to  solve  several  simple  struc¬ 
tures  subjected  to  static  or  dynamic  loads.  It  is  found  that 
even  with  the  use  of  a  relatively  small  number  of  elements  of 
the  type  chosen  in  this  study,  the  results  agree  fairly  well 
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with  solutions  given  by  other  investigators  (whenever  such 
solutions  are  available) .  Results  also  indicate  that  both 
material  and  geometric  nonlinearities  have  a  significant  effect 
in  the  total  deformation  of  structures  which  consist  of  beam 
and  plate  elements.  In  the  case  of  plate  structures  subjected 
to  large  bending  effects,  it  is  found  that  tension  in  the  middle 
surface  of  the  plate  causes  a  significant  stiffening  and  an 
appreciable  change  in  the  vibrational  characteristics  of  the 
structure . 

In  conclusion,  the  feasibility  of  the  numerical  procedure 
developed  in  this  study  has  been  demonstrated  by  the  solution 
of  some  simple  problems.  The  computer  program  generated  as  a 
result  of  this  study  should  be  a  useful  tool  to  structural 
analysts.  However,  further  numerical  experimentation  may  be 
needed  before  the  full  potential  of  the  program  is  realized. 
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Figure  1  -  Typical  Beam  Element 
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Figure  2  -  Plate  Element  with  In-Plane  and  Bending  Actions 
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Figure  3  -  Plate  Element  with  In-Plane  Action  Only 
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Local 


Figure  4  -  Plate  Element  with  Bending  Action  Only 


Figure  5  -  Element  Displacements  Corresponding  to  In-Plane  Action 
(i,  j,  k,  £  represent  the  four  nodes  of  the  plate) 
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Figure  6  -  Element  Displacements  Corresponding  to  Bending  Action 
(i,  j,  k,  £  represent  the  four  nodes  of  the  plate) 


Figure  7  -  General  Stress-Strain  Curve 
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Figure  8  -  Tril inear  Stress-Strain  Curve 
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Figure  9  -  In-Plane  Stressed  Plate  Element  with  8  D.O.F. 
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Figure  10  -  Comparison  of  Central  Deflections  for  a  Clamped 
Square  Plate  with  Concentrated  Load  P  at  Center 
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Figure  11  -  Force-Deformation  Relationships  for  Clamped  and 
Simply  Supported  Square  Plates 
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Figure  12  -  Force-Deformation  Curves  for  First  and  Second  Order 
Static  Analyses  of  a  Plate  Supported  by  Edge  Beams 
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Supported  by  Edge  Beams--Time-History  of  Center  Deflection 

2 

for  Load  Intensities  of  1  and  10  lb/ in 
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Figure  14  -  Second  Order  Elastic  Dynamic  Analysis  of  a  Plate  Supported 
by  Edge  Beams — Time-History  of  Center  Deflection  for 

Load  Intensities  of  1,  2,  5  and  10  lb/in 
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Figure  15  -  First  Order  Elastic  Dynamic  Analysis  of  Simply 
Supported  Square  Plate  with  Lateral  and 
In-Plane  Compressive  Forces — Time- 
History  of  Center  Deflection 
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Figure  20  -  Three-Dimensional  Structure  Consisting  of  In-Plane 

Stressed  Plate  Elements 
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Figure  21  -  Elastic  and  Plastic  Dynamic  Analysis  of  3-D  Structure 
Composed  of  In-Plane  Stressed  Plate  Elements--Time-History 
of  Deflection  in  the  Direction  of  F 
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APPENDIX  A 


STRAIN  AND  STRESS  COMPONENTS  DUE  TO  IN- PLANE  ACTION 


Strain  and  stress  components  to  be  used  in  Eq.  (53)  corre¬ 
sponding  to  the  displacement  functions  defined  by  Eqs.  (51),  and 
(52)  are  presented  in  this  appendix. 

A.l  Strain  Components 

The  strain  components  corresponding  to  the  in-plane 
action  may  be  obtained  by  substitution  of  Eqs.  (51) ,  and  (52) 
into  the  following  equation 


The  results  are  as  follows: 
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A. 2  Stress  Components 

The  corresponding  stress  components  may  be  obtained  by 
substitution  of  the  above  results  into  the  following  equation. 


[0 .  ]  =  a 

1  z 


e  +  vs 
X  z 


vs  +  e 
X  z 


The  results  are  (all  terms  on  the  right-hand  side  of 

2 

the  equations  are  to  be  multiplied  by  the  factor  E/l-v  ) : 


tal]  = 


11  2 
-  i  (2ri  -  3ti  +  1) 

a. 


~  (2n3  -  3n2  +  l) 
a 


g  (1  -  O  ( 6n 2  -  6n) 


c-  (2?3  -  3?2  +  1) 


1  3  9 

g  (20  -  30  +  1) 


£  (l  -  n)  (6^2  -  60 

Cl  _ 


I  (n3  -  2n2  +  n)  -  ~  U3  -  2?2  +  o 


~  (n3  -  2n2  +  n)  -  |  (£3  -  2?2  +  O 
X  E- (l  -  O  (3n2  -  4n  +  l)  +  (l  -  n) 


-45+1)] 
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“V 


i  (2ri3  -  3n2  +  l) 

d 

J  ( 2n 3  -  3n2  +  l) 
g-  S(6n2  -  6n) 


[a9J 


g  (2S3  -  3C2) 

13  2 

£  (2£J  ~  30 

-  I  (1  -  n)  (6?2  -  60 


[a 


[a 


[a 


15 


] 


-  !  (n3  -  2n2  +  n)  -  ^  (?3  -  K2) 

-  (n3  -  2n2  +  n)  -  g  (?3  -  S2) 

A  [-E,  ( 3ri 2  -  4ti  +  1)  +  (1  -  n)  (3£2 

-  |  (2n3  -  3n2) 

v  3  _  2> 

“  g  (2n  -  3r|  ) 
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-  g  n(6^2  -  60 
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[a17! 


|  <n3  -  n2)  +  gi  <;3  -  ?2) 


v~  (n3  -  i2>  +  £  <53  -  52) 


>■  t-C  <3n2  -  2n)  +  n(352  -  25)  ] 


'°19> 


I  (2n3  -  3n2) 


v  nJ  o  2, 
-  (2 n  -  3n  ) 


£  (1  -  5) (6n2  -  6n) 


[a 


21 


£  (2^3  -  3^2  +  1) 

13  2 
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I  n(6?2  -  65) 
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APPENDIX  C 


TYPICAL  INPUT  FORMAT 


C.l  STATIC  ANALYSIS-FIRST  ORDER 

STRUCTURE  ANALYZED  IN  SECTION  6.4  OF  THIS  REPORT 


THIN  PLATE  ELASTICALLY  SUPPORTED  8Y  EOGE  BEAMS-**  PLATE  AND  8  BEAM  ELEMENTS 
STATIC  ANALYSIS 


SPACE 

FRAME  FIRST 

OROER 

number 

OF 

JOINTS  8 

NUMBER 

OF 

NOOE  POINTS  1 

NUMBER 

OF 

BEAM  ELEMENTS  8 

NUMBER 

OF 

PLATE  ELEMENTS  4 

JOINT 

COORDINATES 

1 

0.0 

0.0 

0.0 

2 

15.0 

0.0 

0.0 

3 

30.0 

0.0 

0.0 

4 

0.0 

0.0 

15.0 

5 

15.0 

0.0 

15.0 

6 

30.0 

0.0 

15.0 

7 

0.0 

0.0 

30.0 

8 

15.0 

0.0 

30.0 

NODAL 

COORDINATES 

9 

30.0 

0.0 

30.0 

MEMBER  INCIDENCES 


NUMBER  OF  CONSTRAINEO  VERTICES  9 
1111010 
2101010 
3111010 
<*101010 
5101010 
6101010 
7111010 
8101010 
9111010 

MEMBER  PROPERTIES 

-  1  8  .212000E+01  .00000 OE *00 

0000100009  0.25  1.0 

MOOULUS  OF  ELASTICITY  30.00 
NU  0.30 
LOADING 

1  56.25 

2  112.50 

3  56.25 

<«  112.50 

5  225.00 

6  112.50 

7  56.25 

8  112,50 

9  56.25 
EJECT 


. 440  000E  +  00 


. 420  0  00E-01  1.0 
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C.2  STATIC  ANALYSIS-SECONO  OROER 

STRUCTURE  AnALYZEO  IN  SECTION  6.6  OF  THIS  REPORT 


Thin  plate  elastically  supported  by  eoge  beams-6  plate  and  a  bean  elehents 

STATIC  ANALYSIS 

space  frame  second  orcer 
nuhber  of  joints  8 

NUMBER  OF  NODE  POINTS  1 
NUMBER  OF  BEAM  ELEHENTS  I 
number  of  plate  elements  6 
JOINT  COORDINATES 


1 

0.0 

0.0 

0.0 

2 

15.0 

0.0 

0.0 

3 

30.0 

0.0 

0.0 

6 

0.0 

0.0 

15.0 

5 

15.0 

0.0 

15.0 

6 

30.0 

0.0 

15.0 

7 

0.0 

0.0 

30.0 

8 

15.0 

0.0 

30.0 

nodal 

COORDINATES 

9 

30.0 

0.0 

30.0 

MEMBER 

INCIDENCES 

112 

2  2  3 

3  3  6 

<•6  9 

5  9  6 

6  S  7 

7  6  7 

8  16 

112  5  6 

2  2  3  6  5 

3  6  5  8  7 

6  5  6  9  6 

NUMBER  OF  CONSTRAINED  VERTICES  9 
1111010 
2101010 
3111010 
6101010 
5101010 
6101010 
7111010 
8101010 
9111010 

MEMBER  PROPERTIES 

1  8  .2120Q0E+01  *0000  ODE *00  .660000E+00  .  620000E-01  1.0 


0000100006 

0.25 

MODULUS  OF 

ELASTICITY 

30.00 

NU  0.30 

LOADING 

1 

562.50 

2 

1125.00 

3 

562.50 

6 

1125.00 

5 

2250.00 

6 

1125.00 

7 

562.50 

8 

1125.00 

9 

562.50 

EJECT 
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C3.  ELASTIC  DYNAMIC  ANALYSIS-SECOND  ORDER 

structure  analyzed  in  section  6. a  of  this  report 


THIN  PLATE  ELASTICALLY  SUPPORTED  BY  EDGE  AEAMS-4  PLATE  AND  S  BEAN  ELEMENTS 

ELASTIC  DYNAMIC  ANALYSIS  LTYPE  1 

SPACE  FRAME  SECOND  ORDER 

NUMBER  OF  JOINTS  A 

NUMBER  OF  NODE  POINTS  1 

NUMBER  OF  BEAM  ELEMENTS  A 

NUMBER  OF  PLATE  ELEMENTS  4 

JOINT  COORDINATES 


1 

0.0 

0.0 

0.0 

2 

15.0 

0.0 

0.0 

3 

30.0 

0.0 

.  0.0 

4 

0.0 

0.0 

15.0 

5 

15.0 

0.0 

15.0 

6 

30.0 

0.0 

15.0 

7 

0.0 

0.0 

30.0 

A 

15.0 

0.0 

30.0 

nodal  coordinates 

9 

30.0 

0.0 

30.0 

MEMBER  INCIDENCES 

1 

1  2 

2 

2  3 

3 

3  6 

4 

6  5 

5 

9  A 

6 

A  7 

7 

4  7 

A 

1  4 

1 

1  2 

5 

4 

2 

2  3 

6 

5 

3 

4  5 

A 

7 

4 

5  6 

9 

A 

NUMBER  OF  CONSTRAINED 

VERTICES 

9  • 

1111010 

2101010 

3111010 

4101010 

5101010 

6101010 

7111010 

aioioio 

9111010 

LUMPED  MASSES  AT  VERTEX  POINTS 
1  9  >1010006-02 

MEMBER  PROPERTIES 

1  A  .212000E+01  >0000  006  +  0  0  .440000E  +  00  .420000E-Q1  1.0 


0000100004 

0.25 

MODULUS  OF 

ELASTICITY 

30.00 

NU  0.30 

TIME  PARAMETERS 

0.0 

2.7 

0.0001 

LOADING 

0.0 

2.7 

1 

562.50 

2 

1125.00 

3 

562.50 

4 

1125.00 

5 

2250.00 

6 

1125.00 

7 

562.50 

A 

1125.00 

9 

562.50 

EJECT 
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